Experimental mathematics on the magnetic 
susceptibiUty of the square lattice Ising model 



S. Boukraa^, A. J. Guttmann*, S. Hassani^, I. Jensen^-, 
J.-M. MaillardI , B. Nickel+ and N. Zenine§ 

fLPTHIRM and Departemcnt d'Acronautique, Universite de Blida, Algeria 

I ARC Centre of Excellence for Mathematics and Statistics of Complex Systems 
Department of Mathematics and Statistics, The University of Melbourne, 
Victoria 3010, Australia 

§Centre de Recherche Nucleaire d'Alger, 2 Bd. Frantz Fanon, BP 399, 16000 
Alger, Algeria 

II LPTMC, Universite de Paris, Tour 24, 4eme etage, case 121, 4 Place Jussieu, 
75252 Paris Cedex 05, France 

+ Department of Physics, University of Guelph, Guelph, Ontario NIG 2W1, 
Canada 

E-mail: maillardaiptmc.jussieu.fr, maillardaiptl.jussieu.fr, 
tonygams .unimelb . edu. au, I . Jensenams .unimelb. edu. au, 
njzenineayahoo .com, boukraaamail .univ-blida. dz 

Abstract. We calculate very long low- and high-temperature series for the 
susceptibility x of the square lattice Ising model as well as very long series 
for the five-particle contribution x^^^ ^'^'^ six-particle contribution x'^'- These 
calculations have been made possible by the use of highly optimized polynomial 
time modular algorithms and a total of more than 150000 CPU hours on computer 
clusters. The series for x (low- and high-temperature regime), x'^^ X^^^ ^'''e 
now extended to 2000 terms. In addition, for x'^'i 10000 terms of the series are 
calculated modulo a single prime, and have been used to find the linear ODE 
satisfied by x'^' modulo a prime. 

A diff-Pade analysis of the 2000 terms series for x'^' s-r^i x'^^ confirms to 
a very high degree of confidence previous conjectures about the location and 
strength of the singularities of the n-particle components of the susceptibility, up 
to a small set of "additional" singularities. The exponents at all the singularities 
of the Fuchsian linear ODE of x'^^ and the (as yet unknown) ODE of x'^' a^re 
given: they are all rational numbers. We find the presence of singularities at 
w = l/2 for the linear ODE of x''^', and w)^ = i/g for the ODE of x^*^', which are 
not singularities of the "physical" x'^^ and x'^'i that is to say the series-solutions 
of the ODE's which are analytic at u) = 0. 

Furthermore, analysis of the long series for x'^' (and x'^') combined with the 
corresponding long series for the full susceptibility x yields previously conjectured 
singularities in some x'"' i ^ 7. The exponents at all these singularities are 
also seen to be rational numbers. 

We also present a mechanism of resummation of the logarithmic singularities 
of the x^"^ leading to the known power-law critical behaviour occurring in the 
full Xi and perform a power spectrum analysis giving strong arguments in favor 
of the existence of a natural boundary for the full susceptibility x- 
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1. Introduction 



The magnetic susceptibility of the two-dimensional Ising model has attracted the 
attention of many of the mathematicians and physicists who work (or worked) in the 
area of equilibrium statistical mechanics for more than 60 years, ever since Onsager's 
celebrated solution of the free-energy of the two-dimensional Ising model [1] . 

Let aij be the spin at lattice site The two-point correlation function is 

defined as 

C{M,N) = {ao,o(TM,N), (1) 
and the magnetic susceptibility is given by 

M N 

where the magnetisation M is zero for T > T^. and, in the isotropic case, M = 
(1 - s-'*)i/8 for T < Tc (in the scaling limit M, N oo), where s = sinh(2J//cT) . 

To this day, while we still lack a closed-form solution for the susceptibility, we do 
have a great deal of associated knowledge. While there have been many breakthroughs 
in the study of this problem, we would like to mention five in particular. 

Firstly, in 1976, Wu, McCoy, Tracy and Barouch [2] showed that the susceptibility 
can be expressed as an infinite sum of contributions, known as n-particle contributions, 
so that the high-temperature susceptibility is given by 

kT-XHiw) = J2x^^^+'Hw) = i.(l-s*)i.^x^2„+i)(^) 
and the low-temperature susceptibility is given by 

kT-XLiw) = = (l-l/s^)J-5^x^2")(^i;) (4) 

in terms of the self-dual temperature variable w = ^s/{l + s"^). 

As is now well known [2], the n-particle contributions have an integral 
representation and are given by the {n — l)-dimcnsional integrals [3, 4, 5, 6], 

-. n — 1 «27r II " 2 



where S 



and 



= n = ^^^^((^'-^^v^^-v^, (6) 

-L Xj, 3j 7 



1 < i < J < n 



D(n) _ ^ + nr=l 

1 - ni=i xi 

The Fermionic term G'"^ has several representations [4]. 



(7) 
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with 

Xi = , (8) 

1 - 2wcos(^j) + Y(1 - 2u;cos(^i))^ - 4«;2 

J {1 - 2w cos{(l)i))^ - Aw'^ j=i 

valid for small w and, elsewhere, by analytic continuation. The variable w corresponds 
to small values of s as well as to large values of s. In discussing analytic properties 
of x^"' we will often refer to the principal disc by which we mean |s| < 1 for x*-^""*"^-* 
and |s| > 1 for x^^"''- For simplicity in writing we will however always refer to the 
principal disc as |s| < 1 for the generic x^"^ and leave it to the reader to understand 
the correct high/low temperature distinction. 

Next, in 1996, Guttmann and Enting [7] gave compelling arguments (though not 
a proof) that the anisotropic Ising susceptibility was in a different class of functions to 
that of most functions of exactly solved lattice models. In particular, both the Ising 
free-energy and magnetisation arc known to be holonomic functions (i.e. differentiably 
finite or D- finite functions), while the susceptibility, they argued, was not. This is an 
important clue as to why the susceptibility is, mathematically, a much more difficult 
problem than the free-energy or magnetisation. 

In 1999 and 2000, Nickel [3, 4] suggested that the isotropic Ising susceptibility 
possessed a natural boundary on the unit circle \s\ = 1. While again not providing a 
rigorous proof. Nickel's arguments were most compelling. Note that functions with a 
natural boundary cannot be D-finite. 

Then in 2001, Orrick, Nickel, Guttmann and Perk [8] presented a polynomial 
time algorithm for the generation of the coefficients of the series expansion of the 
susceptibility, in fact in time O(iV^) for a series of N terms. From an algebraic- 
combinatorics viewpoint, a polynomial time algorithm is considered a solution. 
Naturally, a closed form solution is preferable, both for elegance, and so that the 
analytic structure can be readily extracted. Furthermore, Orrick et al. [8] gave a 
detailed asymptotic analysis, based on a series expansion of some 323 terms in both 
high- and low-temperature expansions, and most crucially, short distance correlation 
functions as series in r = (1/s — s)/2, the temperature deviation from the critical 
value. Various exponents and amplitude parameters were extracted to unimaginable 
accuracy, typically 30 or more digits, and these answered most of the outstanding 
questions as to the structure of the scaling functions. 

Orrick et al. [8] also addressed the question of the implications of the unit 
(complex s-plane) circle as a natural boundary [3, 4]. Isolated singularities can 
obviously be "seen" and indeed Orrick et al. needed to subtract or otherwise suppress 
the effect of a number of such singularities, including a contribution from x^^^ to 
achieve the; accuracy they did in their critical point analysis. They also presented a 
heuristic argument and calculation to show that the n ^ oo accumulation of the x^"^ 
singularities implied that expansions in r would only be asymptotic. Their prediction 
of the rate at which expansion coefficients diverge was qualitatively confirmed by 
numerical results on a toy model, subsequently called the diagonal Ising model by 
Boukraa et al. [9]. 

In 2004, Zenine, Boukraa, Hassani and Maillard [10] made an important step 
towards the understanding of the three-particle contribution x^'^-' , based on a novel and 
powerful method for series expansions [11], followed by an inspired computer search 
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for the underlying linear ODE generating the series expansion. They obtained the 

Fuchsian linear ODE for x*-"^-*. In 2005, the same group [12] found the Fuchsian linear 
ODE for x*^**) by similar methods. An important observation coming out of the x^^'' 
work was that there were singularities that had not been predicted by Nickel [3, 4] and 
that the x''"-' were in fact much more complicated functions than had been imagined. 
This gave considerable urgency to finding new results for higher order x^"^- 

Since finding the linear ODE for x*^"^ with n > 4 is clearly a huge task, another 
strategy to get some kind of handle on the analytic structure of x^"^ was considered. 
A set of simplified integrals were introduced [13, 14] by making the integrand of x^""^ 
simpler and/or by restricting the multiple integral to an integral over a single variable. 
The hope was that these integrals will preserve some (or hopefully all) of the singularity 
structure of the full problem. From these model integral investigations, a reasonably 
consistent picture of the singularities emerged and a large set of candidate singularities 
for those actually occurring in the linear ODE's of x''"'' was obtained. 

In this paper we return to the problem of finding linear ODE's for x'^"'' (for 
n > 4) or, failing that, to obtaining long series that can be analysed by other means. 
To save space and repetition, we will drop the important adjective linear before ODE 
in the following, but all mention of ODEs should be taken as a reference to a linear 
Fuchsian ODE, unless otherwise stated. A strong motivation for obtaining such long 
series for the full susceptibility x, and the n-particle contributions x''^'' ^nd x''^'' is 
not to improve already known numerical calculations (critical exponents, amplitudes, 
...), but to shed some light on important, and not yet understood, physical problems 
such as the structure of the ODE's they satisfy, the implication of new singularities for 
the natural boundary, heuristic asymptotics of Orrick et al. [8] , and the mechanism of 
resummation of logarithmic singularities. 

Indeed, in all previous analyses of the full x ^nd the individual x'^"\ ^ point 
appeared that was left unresolved [15]. This issue is the power/log behaviour of each 
X^"^ at the singular points which is not the same as the behaviour of the full x- That 
is to say, the asymptotic behaviour of the full susceptibility x is (with r = |(l/s — s)) 

X ~ ct. |r|-^/^ +ct. +ct.\T\^/^ +0(|T]l0g|T|), (10) 

whereas: 

X(") - ct. |t|-^/4 +ct. |T|i/4(log|r|)"-i +0(log|r|"-2), (11) 

where ct. denotes constants. In order to resolve the issue of how the individual terms 
of the form |r|^/^(log |t|)"~"^ occurring in x^"' combine to give a constant in the full 
susceptibility, we require a better understanding of each x*-"-* rather than getting the 
full X to higher accuracy. With the complete solutions for x^^' and x'^\ long series 
expansions for the higher terms x^"^ and for the full x will allow us to completely 
resolve this issue. 

Our achievements, with regard to series generation, are that we have extended 
the series for the full Xj X^^^ and x^^^ to about 2000 terms§. We also have series 
for x^'''''j modulo a single prime, to 10000 terms, and this has enabled us to find the 
associated linear ODE modulo a prime and get a complete picture of the singularities 
and exponents in that case. Most of the details of this x'^^ linear ODE analysis will 
be given in a future paper. 

In this paper we report in Section 2 on how the generation of the long series 
has been achieved and in Section 3 we describe an efficient method for finding the 

§ More precisely, for x*-^^ we have 3260 coefficients in w or 1616 non-zero terms (that is in a; = w'^). 
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associated linear ODE's and give some results from the linear ODE mod prime analysis 

of x^^y Section 4 is devoted to a floating-point diff-Pade analysis of the 2000 terms 
exact series for x^^^ . There is no such numerical work on long series in the literature 
and so we have had no a priori knowledge about what to expect. In view of this, we 
have been strongly motivated to make a detailed comparison between the numerical 
and exact x^^' results. The comparison has been instructive and we now have a 
reasonably good basis for judging the efficacy and reliability of the diff-Pade method. 
Specifically, there arc clear patterns showing which singularities and exponents can 
be obtained and which will be missed. We report also on the floating-point diff-Pade 
analysis of the 1600 term exact series for x^^^ where only limited exact results are 
known based on our previous [13, 14] toy model investigations. Section 5 is devoted 
to a study of the singularities that occur in the higher x^"'\ n > 7, by performing 
a diff-Pade analysis on the full x series from which the small n contributions X^"\ 
n < 7 terms have been removed. We address in Section 6 the "logarithm summation" 
problem discussed above, and resolve it. A simplified model is introduced to show 
how the subdominant singularities in x^"^ , that individually vanish as the temperature 
approaches its critical value, sum to yield the known finite amplitude at the critical 
point. In Section 7 we describe the fast Fourier transform (FFT) that was used by 
Orrick et al. [8] as a diagnostic tool in the study of series singularities. Here we use it 
to confirm the absence of certain singularities on the complex s-plane unit circle and 
confirm the existence of a natural boundary for x- 

Except for Appendix A on Fuchs' relations, the Appendices revisit and extend 
the Landau singularity analysis [13, 14, 16] of x^"\ In Appendix B we prove that in 
the absence of the "Fermionic factor" {G{nj)'^ in (5), the singularities found in [13, 14] 
are exhaustive. We also prove that none of these singularities, beyond those found 
by Nickel [3], can lie on the principal s-plane unit circle. This dispels any hope that 
there might be singularity cancellation and that, as a result, the unit circle may not 
be a natural boundary. 

The presence of the Fermionic factor in the x^"'' integral (5) is too complicated to 
deal with in its entirety, and Appendices B and D provide a more limited perspective. 
We describe in Appendix D what happens when all, but a single term of the complete 
Fermionic product, is dropped from the integrand of x'^"''- The linear ODE describing 
the series in that case has additional singularities. On the other hand we show that 
there are no such Landau singularities [13, 14] on various analytic continuations^ of the 
-^(") T^lic conclusion seems to be that the singularities of the integral representation 
of x^"'' do not identify with but are just included into singularities of the linear ODE 
associated with a x^"' . In Appendix E we give a power counting argument to determine 
the singularity exponents of one class of Landau singularities. The close analogy of 
these to the principal s-plane unit circle singularities, and the simple formula for the 
exponent values, motivate us to call them the "normal" exponents. Other exponents 
are known to be present based on the exact and diff-Pade analysis in Sections 4 and 
5, and, as a tentative terminology, we denote these other exponents as "anomalous". 

^ The Landau singularity analysis [13, 14] is local. This amounts to keeping no track of which local 
square-root branch we are on for each singularity. 



Extreme series 



6 



2. Extremely long series for the susceptibility 

2.1. The full susceptibility X 

Orrick et al. [8] used an algorithm of complexity 0(A^^) to obtain the first N 
coefficients of the susceptibility series. As a result, they generated and analysed 
series with more than 300 terms in both the high- and low-tcmpcrature regime. This 
remarkable polynomial growth of the algorithm is due to the use of quadratic partial 
difference equations (see [8]) which are finite difference Painleve equations [17, 18]. 

The susceptibility x is obtained via (2) with this quadratic partial difference 
system of equations providing a means to calculate the two-point correlations 
C(m,n) efficiently, and simultaneously, for both high- and low-temperature series. 
A susceptibility scries of N terms requires C{m, n) on the octant m, + n < 2N, 
m < n, deduced as series from the quadratic difference equations, with the diagonal 
C{n,n), n < N, as initial value data. The latter can be obtained from either a 
Toeplitz determinant expression [17] or a set of non- linear recursion relations [19]. 
The calculation of C{m, n) on a particular site requires the multiplication of the 
C{m', n') series of length 0(A'') for the same and/or neighbouring sites. If the series 
multiplication is done as a naive product then the complexity of the susceptibility 
determination is that of 0{N^) integer multiplications on 0(A^^) sites. The Maple 
code used by Orrick et al. [8] ran in a time approximately proportional to N^. This 
was as expected since the integer (digit) size grows linearly in N and naive integer 
multiplication scales as N^. There was no attempt to improve on this code as the 323 
term series obtained was deemed entirely adequate. 

Obvious improvements can be made. First, we know that x when expressed in the 
variable s/2 (high temperature) or l/(4s'^) (low temperature) is a series with integer 
coefficients. This implies that if the series generating program is run using modular 
arithmetic [20] with a number of different primes, then the full series coefficients can be 
reconstructed from the residue series via the Chinese remainder theorem. The number 
of primes necessary is easily estimated. Since |s| = 1 is the singularity boundary for 
both high and low temperatures, the high-temperature series coefficients in s/2 grow 
with the number of terms iV as 2^, while at low temperature the series coefficients in 
l/(4s^) grow as 4^. The latter expressed in numbers of bits is 2N and with primes 
of bit length 15 that we typically use, we arrive at ~ 2N/15 as the required number 
of primes. 

The modular arithmetic approach leaves the complexity of the x series generation 
to TV terms at 0{N*) multiplications for a single prime. For the full calculation 
that requires 0{N) primes the complexity is O(iV^). However, once the generating 
program is running in integer mod prime mode, it is easy to implement series 
multiplication via fast Fourier transform (FFT) [20]. This changes the series 
multiplications from 0{N^) to 0(iVlog(A^)) and gives as our final estimate for the 
X generation complexity 0(A''^ log(A'')). In the next subsection we show how similar 
complexity improvements in the generation of the x*-"' series have been achicived. 

For the extension of the x scries to 2000 terms as reported below, we translated 
the Maple code from Orrick et al. [8] to Fortran^ and ran it in integer mod prime 
mode. We also incorporated the FFT multiplication of series. Finally we changed the 
C(n, n) initialization to use the .limbo and Miwa recursion [19]. 

^ In preliminary studies we obtained 1600 terms of the liigli and low temperature expansions of x 
with highly optimised C-h+ programs. 
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With these changes the calculation on the APAC (Australian Partnership for 
Advanced Computing) SGI Altrix cluster with 1.6GHz Itanium2 processors using 280 
primes took about 240 CPU hours in total (we note that the algorithm without FFT 
multiplication would require about twice the above amount of CPU time). 

The resulting series can be found on the wcb-pagc [21] where we give the first 
2000 coefficients of the high- and low-temperature expansions in the variables u = s/2 
and V = l/(4s^), respectively. In addition we also give the expansion for xh and xl in 
the self-dual variable w ~ is/(l + s^). The scries in the u or v variables, or in the w 
variable, are probably the best as far as computer encoding and modular calculations 
are concerned since all the coefficients are integers rather than rational numbers. For 
example, for any n > 2, x'"^ = 2"w" [1 -|- An'^w'^ + 0{w^)]. 

2.2. The contributions x^^^ o-'^d x^^^ 

The previous longest series available for x^^^ and x^^^ can be found in [8] where the 
first 182 and 140 terms respectively of the series in s are listed. We have extended 
these series dramatically - two extreme examples being the full integer series for x^^^ 
to 2000 terms and the series for x'^'^' modulo a single prime to 10000 terms. The 
latter extension is only possible because the complexity order for a mod prime series 
of length N has been reduced to 0(A'^^ log(A'')). We outline in this section how this 
reduction has been achieved. 

Our method for evaluating the integral expression (5) for x^"'' remains as 
described in earlier publications [11, 12]. We first convert (5) back to an n-fold 
integration with the explicit phase constraint 27r5(^ 0i) now in the integrand. This 
constraint is then replaced by the equivalent exp(ifc '^(pi)^ thus decouphng all 
integrations at the expense of a sum over the Fourier integer^ k, where the sum 
extends from /c = — oo to oo. 

Next we expand all denominator factors in the integrand of x^"''; thereby 
converting the integrand into a sum of n-fold products J| j/j x"* . Each 0i integration 
then picks out the k*^ Fourier coefficient of j/ix"', namely w''^' +"<+i . a{k,ni), where 
a{k,ni) is proportional to a hypergeometric function in the variable 16w^. 
The net result of all these operations is that we have replaced the x^"^ continuum 
integration by a nested summation of products of hypergeometric functions. 

The complexity of this calculation is of some order that we can now easily 
determine. If we want a series of length A'", then the Fourier k sum can be restricted 
to 0(A'') as can all the hypergeometric function series. The evaluation of products 
of series of length 0{N) is either an 0(A^^) calculation if done as a naive product 
or 0{N\og{N)) if done by FFT. These two operations, namely the k sum and the 
series multiplication, are inherent to our method and give an irreducible minimum 
complexity of either 0{N^) or 0{N'^log{N)). The only place where we have some 
freedom to reduce the complexity of our x^") series evaluation is in the number of 
summations that are required for the expansion of the denominator factors in the 
original x^") integrand. The total number of products Yl Ui x^' must be no more than 
0{N'^) to keep the overall complexity at 0(A^'' log(A')) and this, in turn, implies the 
denominator expansion must be limited to two independent summations. 

^ Trigonometric functions in the integrand such as sin((^, — (/)j)/2) are dealt with by expanding them 
into a sum of phases exp(i^mj(/ii) with each rm some small integer. These can be incorporated 
into shifts ktl>i ^ {k + mj) • 4>i and do not change the form of the integrand. Since such shifts also 
do not change the complexity order of the x^"' calculation they will not be considered further. 
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This limitation on the summations immediately shows that the product form 

(6) for the Fcrmionic factor G*-"-* is not appropriate and alternatives must be used. 
The useful formulae for high- temperature series are those given in [3], in particular 
equation (5) for some low order and equation (10) for the general case. For 

low-temperature series we note that G*^^"^ has been shown to be a Pfaffian [22] in 
the hij defined in (6). By a rearrangement of terms in the Pfaffian one finds that the 
analog of equation (10) in [3] is 

o-(2n) _ (G^'^"^)^ 1 , , 

^ - -{2njr - (2;i)'^i2/i2„-i,2„ (12) 

n — 1 ^ 
X I I [h2m,2m+l h2m-l,2m+2 +77; ^h2m-l,2m h2m+l,2m+2 

V (2m) 

m=l 

Use of label interchange symmetry, hij = —hji, allows one to combine terms in (12) 
further. For the first few low order terms we have 

H^^^ = -^(/Il2/l2l), 

-H"^^^ = ^ {hl2 h2l){hs4 /I43) - ^ C^lS /l23 h34h4l), 

H^^^ = --^(/Il2/i2l)(/l34/l43)(/l56/l65) + ^ (/ll2 /i23 /i34 /l4l)(/l56 /I65) 

- 7 {hi2 h23 /i34 h^a hei), (13) 
6 

written in an obvious cyclic form. It is these expressions, in particular the one for 
H'^^^ , that are the starting point of our discussion of the reduction in the denominator 

expansion summations. The most complicated term contributing to x*-^-* has an 
integrand that contains the last term of H'-^^ in (13) and thus the seven denominator 
factors 

(1 - Xi X2 X3 X4 X5 Xe)~^ (1 - Xi X2)~^ (1 - X2 0:3)"^ (14) 

X (1 - a;3 0:4)"^ (1 - X4 a;5)~^(l - X5 a;6)~^(l - xe xi)~^ 

Naive expansion of these denominators results in a 7-fold sum clearly a 5-fold excess 
that must be eliminated. The first step towards this elimination is the use of partial 
fraction rearrangement, a trick that was already used in the evaluation of x'^^ and 
X^*) in [11, 12]. The version we use here is based on the identity 

(1- A)-i(l-B)-i(l-G)-i = 

(1 - ABC)-^ [1 - (1 - A)-^ - (1 - B)-i - (1 - G)-i + (15) 
(1 - A)-i (1 ~ B)-i + (1 - B)-i (1 - G)-i + (1 - G)-i(l - A)-^], 

in which we first set A,B,C = X1X2, X3X4, x^xq and then A,B^C = 
X2X3, X4X5, xqXi. The product of these two forms of (15) enables us to replace 
(14) by 

(1 - Xi X2 X3 X4 X5 Xq)~^ (1 - Xi X2)~^ (1 - X2 X3)~^ 

X (1 — X3 x^)'^ (1 — X4 x^)^^ (16) 

plus terms with similar, or fewer, denominators. That is, (16) is now the most 
complicated set of denominators in the integrand contributing to x^^^ With the 
replacement of the first factor by (1 — xi X2 x^ X4 0:5)"^, (16) becomes the equivalent 
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most complicated term in the evaluation of x^^^ . Thus all our subsequent remarks 
apply equally to both x*^^^ and x^^^ . 

Expansion of (16) results in the formal 5-fold summation 

C _ \ ^ m+ni m+ni+p m+p+q m+q+ns m+ns m /lyN 

O — / ^ ^1 .^2 •''3 -^4 -^5 •''6 ' ) 

which is of complexity 0{N^). However the rii and summations can be done 
independently and thus (17) is in fact only of complexity 0{N'^ ■ {N + N)) = 0(A^'*). 
This is not yet an adequate reduction and we can do better by eliminating the ni and 
ns summations entirely via the use of recursion relations. For example, we define the 
111 sum of the xi , X2 pair as 

Ai2{m,p) = xT^"^x^+-^+^ = (18) 

ni=0 r=m 

and note that because m appears only as a limit on a dummy variable sum, the A12 
satisfies the recursion 

Ai2{m,p) = SAi2{m,p) + Ai2{m+l,p), (19) 



dAi2{m,p) = x" 



m m+p 



X 



2 



Only the lower limit on the summation in (18) has been given explicitly. There is also 
an upper limit that depends on the length N of the series in x^^' or x'^' want 
to obtain. Thus Ai2{m,p) vanishes for m large enough, and (19) shows that, as m 
is decreased from this upper limit, each determination of Ai2{m,p) requires only the 
evaluation of a single product and its accumulation into a previously stored result|. 
If we now take it as given that the m summation in (17) is performed in decreasing 
sequence we obtain 

S = 5] (m, p) • x^+P^" ■ ^54 {m,q)-x^, (20) 

which is of complexity 0{N^). 

It remains to be shown that the idea of recursion can be applied onc;c; more, 
ultimately reducing the complexity of the calculation of S' to 0{N'^). For this purpose 
define 

Si23(m,g) = 5^^i2(m,p)a;™+P+«, (21) 

which, for q > 0, can be put in the recursive form 

Bi23(to, q)= Y1 i^^M^^P) + ^i2("^ + xT^"^" 

p=0 

= dBi23{m, q) + Bi23{m +l,q- 1). (22) 
For SBi23{'m,q), which is the first sum in (22), we have 

5Bi23 (m, g) = ^xT x'^+P = xT ■ ^23 (m, 9) , (23) 

p=0 

X Our discussion here is schematic. It is to be understood that the ij>i integrations have been carried 
out and the "single product" being referred to is the product of the two hypcrgeometric function series 
associated with x^^x^'^^. An additional implication is that Ai2{ni,p) has absolutely no functional 
dependence on its subscripts and could equally well be denoted A{m,p). The subscripts have only 
been included to indicate a connection to a particular fa<;tor in the x*-^^ or x'®^ integrand. 
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a single product analogous to 5Ai2 in (19). The new feature in the recursion (22) 
for -Bi23, relative to (19) for A12, is that wc must maintain in storage an entire array 
of elements indexed by q. Furthermore, for each m we must supply, by a separate 
calculation, the q = term 

i?i23(m,0) = ^Ai2im,p)-xT+P. (24) 

Although this does require a sum, the fact that g = is fixed means the contribution 
of this evaluation to the complexity order of S is still only 0{N'^). With the general 
-B123 now given either by (22) or (24), we obtain 

S = Yl BMm, q) ■ A^iim, q) ■ x"^, (25) 

which is the two-fold sum, and thus the 0{N^) result, we were looking for. Note that 
memory requirements are also quite minimal. We need to store the B array which has 
0{N) elements each of which is a series of length 0{N). Thus memory requirements 
also scale as N"^. 

Our initial Fortran coding for the x^^^ and x^^'^ series generation used only naive 
series multiplication and thus was 0{N^) for a single prime. With these programs, 
series for x^^^ to 2000 terms and x^^^ to 3260 terms were generated in about 100000 
CPU hours running 160 primes in parallel. The calculations were carried out on the 
afore-mentioned APAC cluster. Series for x'-^^ to 6000 terms, modulo the single prime 
Po = 32749, were obtained in about 40000 CPU hours using 32 processors on one 
of the VPAC (Victorian Partnership for Advanced Computing) facilities which is a 
Linux cluster based on Xeon 2.8Ghz CPUs. 

When it was observed that the 6000 terms were not sufficient to obtain the linear 
ODE for x^^\ the Fortran codes were modified to include FFT series multiplication 
and series modulo po for x*-^-* to 10000 terms were obtained. 

The total CPU time for x^^^ to 10000 terms was about 17000 hours (6000 terms 
take around 2000 hours using the FFT algorithm on the APAC). 

The calculation was done in parallel using 128 processors. The algorithm is 
straightforward to parallelise because the calculations in the outer most loop, that 
is the sum over Fourier mode integer k, can be done independently for each value 
of k. The only issue is that the time required decreases with k. In order to use 
approximately the same time per processor we simply assign calculations with a given 
k to processors in an alternating pattern such that processor does k = 0, k = 255, 
k = 256, k = 512, k = 513,... while processor 1 does k = 1, k = 254, k ~ 257, 
k = 511, k = 514, . . . and so on up to processor 127 which does k = 127, k = 128, 

k = 383, k = 384, k = 639, This simple assignment scheme ensures a good 

balance with the total time used by various processors differing by less than 4%. 

It is this 10000 terms series that has enabled us to obtain the exact ODE modulo 
Po for x^^\ thus making possible the various comparisons found elsewhere in this 
paper. 

We conclude with two observations on the technical aspects of the mod prime 
calculations. Firstly, a very frequent operation in our codes is the accumulation of 
two products, i.e. a a + be + de. If the variables arc integers m,odulo p then 
the accumulation can be done as a standard 32 bit integer operation with the single 
(Fortran) call a = mod(a + bc + de,p) without overflow provided p < 2^^. It is this 
feature that dictates our choice of primes and in particular po = — 19. Secondly, 
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we follow a recommended practice [23] of loading the floating point FFT routines with 
mod prime variables in "balanced" form. That is, if any mod p variable v is greater 
than p/2 it is loaded as v — p. This eliminates most of the "dc" part of the input 
signal and typically increases the safety margin in the output rounding of float to 
integer by several bits. We keep track of the differences in these rounding operations 
so as to guarantee our programs generate all integers correctly. 

3. Fuchsian ODEs for long series modulo a prime 

3.1. The linear ODE for x'^^^ 

We begin with a remark that applies to all the subsequent discussion, namely that 
there is no single unique ODE that describes any given series. There is a minimum 
order linear ODE that is unique but this typically contains a very large number of 
apparent singularities and can only be determined from a corresponding larger number 
of series coefiicients. In our quest for linear ODE's corresponding to given very long 
series expansions, we are interested in the Fuchsian linear ODE requiring the minimum 
number of coefficients from the series in order to be obtained. In general there are any 
number of intermediate ODE's but unless otherwise required by the context, we will 
call all of these the underlying linear ODE, without distinction. 

The 2000 terms generated for x^^^ and X^^^ are not sufiicient to find the exact 
underlying linear ODE. However an alternative approach is to use mod prime series to 
find the linear ODE modulo a prime. From such a linear ODE we can get singularity 
positions modulo a prime, indicial equations and hence singularity exponents modulo 
a prime, and indeed practically everything that could be obtained from the exact 
ODE but restricted to prime residues. Furthermore if, for example, the modulo prime 
factorization of the head polynomial of the linear ODE yields factors with small 
coefficients, then one can, with almost perfect certainty, conclude that one has all 
the exact singularity locations. Similar remarks would also apply, say, to singularity 
exponentsf. And indeed, because we have been able to find the mod prime linear ODE 
for x^^"^ and because the singularity locations and exponents appear to be simple, we 
are confident that what we report are in fact the exact values. 

Because series generation is expensive we want to be sure our algorithm for 
deducing the underlying ODE requires the fewest number of terms. We report in 
this section on a method [24] that is slightly different from that used for, say, x'^^^ but 
that appears to have a number of appealing advantages. We will report on some of 
our results on the analysis of x^^\ in particular those that affect the number of series 
terms required to find the ODE. Most of the details, such as the factorization of the 
linear ODE, will be left for a future publication. 

An essential constraint on the linear ODE underlying a series S{x) of the type we 
are considering here is that it must be Fuchsian. Specifically this means that a; = 
and X = oo are regular singular points. A form for the linear differential operator 
that automatically satisfies this (0, oo) regularity constraint is: 




M 



D 



m=0 d=0 



f Of course, if there is an ambiguity one can always use additional primes and resolve the ambiguity 
by the Chinese remainder theorem. 
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The omo 7^ condition is the obvious constraint to make x = a regular singular 

point and it is the use of the operator x rather than d/dx that makes analysis 
around x = oo simple. The change of varialsle x = 1/y turns (26) into 

M D ^ 
Lmd = Y.H'^^'i- y"""^ ■ ^)"' "^^o ¥= 0, AMD ¥= 0, (27) 

m=0 d=0 y 

where one can see that the condition for x = oo {y = 0) to be a regular singular point 
is ttMD ^ 0. A simple rearrangement of terms shows that Lmd can also be written 

M D , 

Lmd = T.12^rn,-x''+"^-{-r, (28) 

m=0 d=0 

with the b coefficients being linear combinations of the a. There is no single b coefficient 
analog of Umd 0. We will use (26) exclusively, particularly for the purpose of 
determining the Umd, but this does not preclude transforming to (28) if required. 

To determine the Umd in (26) we demand Lmd{S{x)) = and this yields a set 
of linear equations that we arrange in some well defined order. There exists a non- 
trivial solution if the Nmd x Nmd determinant (with Nmd = {M + 1) • (Z? + 1)) 
corresponding to the chosen ordering vanishes. We test this by standard Gaussian 
elimination, creating an upper triangular matrix U in the process. If we find 
U{N,N) = for some N, a non-trivial solution exists. If A'' < Nmd we set to 
zero all amd in the ordered list beyond N. Of the remaining Umd we set omo = 1, 
thus guaranteeing that x = is a regular singular point, and determine the rest by 
back substitution. We guarantee that x = oo is a regular singular point by choosing 
the initial ordering of elements such that no matter what N is, the first elements set 
to zero will be those from row aod, then row aid, etc. In this way the element Umd 
will never vanish unless all UmD, m < M vanish or the series S{x) has x = oo as an 
irregular singular point. We have not systematically investigated what happens when 
we change the element ordering within these constraints. 

The N for which U (N, N) = is the minimum number of coefficients needed 
to find the linear ODE within the constraint of a given M and D. Obviously, 
N < Nmd = {M + 1) ■ {D + 1). Henceforth, D will always refer to the minimum 
D for which a solution is found for the given M. Then we can define a unique non- 
negative deviation A by TV = Nmd — A. Examples of such constants are given in 
Table 1 based on our analysis of x^^^ and the combination 2 • x^^^ — x^^^- 

A very striking empirical observation arises from Table 1 and has been checked in 
many cases as summarised in Table 2. For reasons we do not understand, there exists 
the linear relationship 

N = A- M + B ■ D - C = (M + 1) ■ (D + 1) - A (29) 

where A, B and C are constants depending on the particular series S{x). For x'^' 
they are A = 72, B = 33, C = 900, while for the combination 2 • x*^^ - X^^^ they 
arc A = 68, B = 30, C = 744 as can be verified from Table 1. Note that (29) has 
no (positive) solution for D if M < 5. Thus B = Mq is the minimum order possible 
for the linear differential operator that annihilates S'(a;)§. Similarly, A = Dq is the 
minimum possible degree and thus we can rewrite (29) in the more definitive form 

N = Do-M+Mo-D-C = (M + 1) • (£) + 1) - A. (30) 

§ Generically B = Mq. However, if the ODE is such that aod = for all d then the constant 
B = Mq — 1, that is, the minimum order minus 1. 



Extreme series 



13 



Table 1. M is the order of the linear ODE, D is the degree of each polynomial 
multiplying each derivative, Nmd = + + 1), N is the actual number of 

terms predicted by (29) as necessary to find an ODE of the given order M , and A 
is the difference Nmd — N. The first five columns gives this data for x'^^ while 
the next five columns gives this data for 2 • x^^) _ ^(3) ^ clearly showing the saving 
in the number of terms needed to identify the linear ODE. 



Terms needed to find x 


Terms needed to find 2 • x^^^ - x^^' 


M 


D 


Nmd 


N 


A 


M 


D 


Nmd 


N 


A 


52 


141 


7526 


7497 


29 


48 


131 


6468 


6450 


18 


53 


137 


7452 


7437 


15 


49 


128 


6450 


6428 


22 


54 


134 


7425 


7410 


15 


50 


125 


6426 


6406 


20 


55 


132 


7448 


7416 


32 


51 


123 


6448 


6414 


34 


56 


129 


7410 


7389 


21 


52 


120 


6413 


6392 


21 


57 


127 


7424 


7395 


29 


53 


118 


6426 


6400 


26 


58 


125 


7434 


7401 


33 


54 


116 


6435 


6408 


27 


59 


123 


7440 


7407 


33 


55 


114 


6440 


6416 


24 


60 


121 


7442 


7413 


29 


56 


112 


6441 


6424 


17 


61 


119 


7440 


7419 


21 


57 


111 


6496 


6462 


34 



The minimum order Mq and degree Dq can also be inferred directly from the ODE 
independently of (30). 

The head polynomial Y^^=o ^'nd ™ (26) can be factored modulo a prime and 
the greatest common divisor of these, from several different Lmdi is the polynomial 
P whose zeros are the "true singularities" of the linear ODE. In all cases we have 
tested, the degree of this head polynomial factor is the Dq in (30). For x'"'''("') 
factor, P{w) = P]j^{w) = ^72(111), is sufficiently simple that we are confident that it 
equals what one would obtain from the exact (not modulo a prime) ODE. We report it 
(cf. (34)) in Section 3.3 in the context of a more general discussion of singularities and 
make extensive use of it in Section 4 for comparison purposes with results obtained 
from floating-point diff-Pade analysis of the 2000 term exact x'^' series. 

The multiplicity of any zero x = Xg of Pdo{x) is the number of linearly 
independent singular functions in the neighbourhood of Xg\^. If wc add to this list 
of multiplicities the number of independent singular functions at a; = and at 
y = 1/x = 0, then the maximum multiplicity is the minimum order of the ODE. We 
have again verified that this agrees with Mg in (30), although there is a subtlety 
to determining the true multiplicities at a; = and y = 0. By our definition 
of Lmd in (26) the indicial equations at these points are, X]m=o ■ P"^ = 0, 
and J2m=o ■ (— p)™ = 0, respectively, and thus of degree M. If there are 
positive integer roots then, even if we factor these indicial equations modulo a prime 
for several Lmd and take the greatest common divisor, we can only be sure that the 
true multiplicity is less than or equal to the degree of the greatest common divisor. 
The ambiguity at a: = 00 is easily resolved by analysing the transformed series S{z) 
where x = z/(l — z/a), thus mapping x = 00 to 2 = a. In all our examples the 
situation at a; = has been unambiguous but we do not see why this would be the 

tl There is no difficulty in finding the roots of the associated indicial equation mod prime. In the 
case of x'^' ™s are again confident that we have all the exact exponents and these are reported in 
Section 4. 
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case in general. 

To interpret the constant C in (30) set M = and define A = Aq in this 
case. Then D = {Dq — 1) Mq — C + Aq — 1. On the other hand we know that when 
M = Mo, the head polynomial of the linear ODE factorizes into two polynomials, 
one of degree Dq giving the true singularities of the Unear ODE and one of degree 
Dapp whose zeros are all apparent singularities. Thus we can write D = Dq + Dapp 
which, combined with the solution D = {Dq — 1)Mq — C + Aq — 1, yields 

C = (Dq-I)- (Mo - 1) - D^pp + Ao - 2, (31) 

giving a direct connection between C and the apparent polynomial of the minimum 

order linear ODE. We remark further that Dapp is related to a (true) singularity 
exponent sum-rule. In the case of x''^'' the exponents we have determined yield 
Dapp = 1384 and hence C = (72-l)(33-l) -1384+Ao - 2 = 886 +Ao > 886, 
which is consistent with the observed C = 900. The details of this sum-rule 
calculation can be found in Appendix A. 

The deviations A observed in Table 1 are quite small and a reasonable 
approximation to (30) is obtained by setting A = 0. Then (30) is both a specification 
of N, the minimum number of series terms needed to get a linear ODE, and a 
relationship between M and D. We can use the latter to eliminate, say, D and 
the former to find the minimum possible A*". The result is that one should be exploring 
the region around (M + 1)/(_D + 1) = Mq/Dq to obtain the minimum N. 

In earlier work [11, 12] we observed a "Russian-doll" structure for the linear 
differential equations for x'^^ and x*^^^ and a similar inherited Russian-doll structure 
on the n-particle contributions of the "diagonal susceptibihty" [9] . We conjecture for 
arbitrary x*-"-* a "strong" Russian-doll structure for the linear differential operator for 
X^"\ which is to say that the linear differential operator for x''"' right-divides the 
linear differential operator for x'"^^'- We can now verify this conjecture on x''^' a-nd 
the results are shown in Table 2. 

Further, a stronger property amounts to saying that we actually have, in the 
decomposition of the linear differential operator for x'"^^\ the linear differential 
operator for x^") occuring as part of a direct sum. Such a reduction was found [10, 
11, 12] for the combinations 6 • x^""''^'' — n ■ x^"'', n = 1 or 2, and we now verify 
this conjecture for the case n = 3. Detailed results that we referred to earlier are in 
Table 1. A summary of all the observed operator reductions appears in Table 2. 

Some operator reduction data in Table 2 calls for explanation. Since the 
differential operator Ly for x^''' [H], acting on the x^'"''' series gives an order 26 ODE 
(see sixth line in Table 2) and the minimal ODE for 6 of order 30, i.e. 

less than 33 (see ninth line in Table 2), one can conclude that this order 30 differential 
operator contains an order 4 differential operator occurring in the known L7. From 
this order 4 differential operator, we focus here on the differential operators of order 
1, which because they have no apparent singularities, are most effective in rediicing 
the number of series coefficients that need to be generated. There are two such order 
1 operators)); their solutions are 

Si = w/{l-iw), S2 = w'^ /{{I - iw)^/l - 16w^), (32) 
We have also found an order 1 operator whose solution 

S3 = wV(l-4w)2 (33) 



tl They already occurred [10] as solutions of L7, the differential operator for x'^^. 
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Table 2. Summary of results for various series. The equation for the size of the 
zero determinant is obtained from fits to data such as that shown in Table 1 for 
X^®). The last five columns are the data for the case that the zero determinant 
size is a minimum. A reduced difference in braces in column 7 appears in those 
cases where a constant is a solution of the ODE and the matrix being tested could 
have been taJsen as size O x {D + 1) rather than (O + 1) x (D + 1). 

The series are the model integrals [14] (see (35)). 



Series 


N = DqM + MqD -C 


M 


D 


Nmd 


N 


A 




IM + ID +1 


1 


1 


4 


3 


1 




IM + 2£) +2 


2 


1 


6 


6 







12M + 7D -40 


11 


17 


216 


211 


5 




7M + - 36 


15 


9 


160 


150 


10(0) 




72 M + 33 £1 - 900 


56 


129 


7410 


7389 


21 


Lrix^'^) 

6;^(3) — 


60M + 26£) +611 


54 


131 


7260 


7257 


3 


12M + 6D - 28 


10 


17 


198 


194 


4 


6^(4) _ 2x^2) 


6M + 7L» - 17 


13 


8 


126 


117 


9(0) 


6x(5) - 3x(^) 


68 M + 30 £> - 744 


52 


120 


6413 


6392 


21 




65M + 28D - 526 


50 


117 


6018 


6000 


18 




64M + 27 £1 -409 


49 


117 


5900 


5886 


14 




lOM + 5£) - 21 


8 


13 


126 


124 


2 




5M + 6£) - 12 


9 


6 


70 


69 


1 




45 M + 17 D -277 


28 


80 


2349 


2343 


6 




26M + 271? -342 


48 


39 


1960 


1959 


1 




145 M + 49 - 1943 


92 


257 


23994 


23990 


4 



is a solution of the order 30 differential operator for 6 x^^^ — 3 x^^^ , but not a solution 
of the differential operator for x^^^- Let us introduce the second order differential 
operator| L2 which simtiltancously annihilates 6*1 and 5*2 and the order 3 differential 
operator L3 which annilates 5*3 as well. When these act on the series 6 x^^^ — 3x^^^ the 
reductions in N are shown in Table 2 (tenth and eleventh lines, labelled respectively as 
L2 and L3). The essential observation is that, although we found in our first computer 
runs that 6000 terms were not enough to get the linear ODE mod po for X^^\ this 
is now more than adequate to get the linear ODE modulo p for x^^^ for as many new 
primes p as we wish to investigate. 

3.2. On the Imear ODE for x^^^ 

We do not have corresponding results for x^^^ and it seems likely that obtaining long 
enough series in this case is beyond our presently available computing resources. We 
base this on the following very crude correspondence: determining the exact linear 
ODE for x*-'^'' required a minimal series of about 220 terms. The exact ODE for x'-^'* 
required about 170 terms (in x = w^). The ratio iV^(4) /A^^o) ~ 170/220 ~ 0.77. 

There is a similar ratio for integrals without the Fermionic factor introduced 

in [14]. From Table 2, we have A^<[,(4) /A^^c) ~ 80/130 ~ 0.62 and A^$(6)/A^$(5) ~ 

X The differential operator L2 was given as 0\ ■ N\ = T\ ■ L\ in eq.(7) of [15]. It is thus the direct 
sum of the order 1 operators associated with 5i and 52 . 
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2000/2400 ~ 0.83. A reasonable guess then might be iV^(6) /iV^cs) in the range 0.9 
to 1.1. Although we generated 10000 terms for x^^^ only about 7500 were actually 
required. Our guess is that about the same number will be required for x'^'. Now our 
codes for x^^^ x^^^ s-^e such that 7500 terms for x^^^ is roughly the equivalent in 
time to 15000 terms for x^^^ and this moans a nmning time for x'^^ about 1.5^ ~ 5 
times in excess of that taken for the 10000 terms of x^^^. This might be reduced 
somewhat if we rely on some direct sum assumptions about the linear ODE for x^^^ 
similar to those we have found work for x^^''- On the other hand, a guess for a series 
length A'' that is too small might leave us with no results whatsoever and thus we 
would probably want to err on the conservative side and require a calculation of x^^^ 
with a run time cost as much as 10 times that for x^^^ . 

3.3. Singularities of the linear ODE for x^^^ and Landau singularities 

From the linear ODE for x*^^^ obtained modulo a prime, one can easily reconstruct the 
singularity polynomials of the ODE as they appear at the highest derivative. These 
polynomials read 

. (1 - 4w)22(l + 4w)iS(l - wf{l + 2wf{l + iw + Aw^f 
(1 + w){l - 3w7 + w2)(l + 2w- ■iu?){l - w - 3w;2 + 4^^) 
(1 + 8w + 20w^ + 15w^ + Aw^){l -7w + - Aw^) 
(l + 4u; + 8w2)(l-2«;). (34) 

All these singularities, except (1 — 2w), have been predicted by the model integrals 

wc introduced in [14]. The models considered in [13, 14] arc integral representations 
(one-dimensional and multidimensional) which belong to the "Ising class" [25]. These 
integrals are holonomic and we obtained the linear ODE's of these sets of integrals 
through series expansions [13, 14]. In [14] a detailed analysis of the multiple integrals 

was performed. These n-fold integrals correspond to removing the Fermionic 
factor (G("))2 in (5), so that 

We obtained (after eliminating the apparent singularities) the following 

(n) 

polynomial factors for the head polynomial of the linear ODE's satisfied by the , 
expressed in terms of Chebyshev polynomials of the first and second kind [14]: 

T2p, (l/2w; + l) = T„_2pi-2p, (l/2u;-l), (36) 

< pi < [n/2], < P2 < [n/2] - pi, 

and the polynomial arising from the elimination of z in: 

"M2«> ^'\2w l-Awz 



/ 1 

l(z)-C/„,_i(— -- 
V 2w 1 



4w — z 



Aw z > 
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ni = pi, n2 = n - pi - 2p2, 

< pi < n, < P2 < [(n-pi)/2]. 

Our motivation was to obtain "good candidates" for the factors in the head 
polynomial of the linear ODE of X*^"\ expecting that the Fermionic factor (G*^"-')^ 
may not introduce "too many additional singularities" . 

There is also the possibility of singularity cancellation. Indeed, in the even simpler 
integral in which the R^""^ factor (7) is dropped from (35), all the singularities predicted 
in [16] are seen; these are the singularities (36) but with the even 2pi replaced by an 
integer that can also be odd. Thus the inclusion of R^"'^ has eliminated a whole class 
of singularities. The presence of (G'"^) , at least in x^^\ has been much less dramatic 
and led only to the one extra factor (1 — 2w) in (34). 

To get a better understanding of the effect of (G'"^) we revisit the Landau 
singularity approach in Appendices B through E but performing the calculations on 
the (n— l)-fold integrals (5) rather than the original (2n — 2)-fold integrals [2] that was 
the basis of calculations in [14] . The new approach detailed in Appendix B confirms our 
previous Landau singularity calculations [13, 14], the difference being in the number 
of integration variables considered. The original representation for x^"^ given in [2] is 
an integral over two sets of (n — 1) independent phases ^j, (j with an integrand that 
is symmetric under the interchange of these sets. Integrating out one set to arrive 
at (5) has obviously broken this symmetry but a vestige of it remains, such that for 
every combination of the (j, (pi leading to a singularity there is another set (pj, d 
obtained by C <-^- (/> interchange that leads to the same singularity. In addition another 
symmetry arises such that, for a given (pi, Q combination, the Landau conditions 
allow the reversed combination —(pi, — Ci explaining the n ^ n — 2m replacement 
singularities seen to occur in the analysis of [14]. Our analysis proves that for the 
4>^-' there are no singularities other than those given by (36) and (37). We also 
show that (G*^"))^ neither reintroduces the singularities cancelled by the i?*^") nor 
leads to further cancellation and thus the singularities of are included in the 

singularities of x^"^. The (G'"^) factor has a dramatic effect on the exponents of 
the singularities and a power counting argument is given in Appendix E to predict 
the exponents at all singularities. Those calculations are not intended to cover all 
possible contingencies and while many predicted exponents are observed there remain 
a number of "anomalous" cases. 

Our analysis confirms the fact that the Landau singularities of x''"\ at least those 
given by are included in the set of singularities of x^"'''^'^^ A; > 1. This is exactly 
what we have found for the linear ODE for x''^''- The first line in (34) corresponds to 
the singularities occurring in the linear ODE for x'''^^- This is a necessary condition for 
X^^^ to be embedded in x^^^. However this comes with an important caveat. We prove 
in Appendix B that none of the embedded singularities can lie on the |.s| = 1 boundary 
of the principal disc of the x'^"'' function defined by the integral (5). We also show 
that of all the different classes of singularities, only those we call Case 2 irreducible 
singularities in Appendix B lie on |s| = 1 of the principal disc. These correspond 
to the singularities derived in [3, 4] and are elsewhere called the Nickelian or circle 
singularities. The importance of this lies in the fact that we cannot expect cancellation 
between different singularities on the principal disc and thus the elimination of |s| =1 
as a natural boundary. 

While some aspects of the effect of (G^"^) have been determined, a complete 
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Landau singularity analysis of x^"^ has not been done; in particular we cannot 
definitively state whether the (1 — 2w) factor in (34) does or does not identify with 
a Landau singularity of x^^^ . However, a toy integral intermediate between and 
discussed in Appendix D and provides an example in which it can be shown 
that the ODE describing an integral has more singularities than the integral. The 
analogy to x^"^ is sufficiently close that we believe it is likely that (1 — 2w) = is not 
a singularity of x^^^ . Clearly it is of interest to know whether there are such additional 
singularities in the ODE for x*-"^ for larger n and whether they are also on \s\ = 1. 

Appendix C lists, for the first few n of the singularities corresponding to 

(36,37) and also derived in Appendix B. In the next section we will see whether these 
singularities appear in the analysis of our long series with exact coefficients. 

4. Diff-Pade analysis of the long series for x^^'' and x*^^^ 

We present, in this section, a diff-Pade analysis^ of our long 2000 coefficient scries 
(these are the actual series coefficients, not the coefficients modulo a prime, for which 
we have longer series) in floating point form, in order to obtain the singularities 
(together with their associated exponents) that should occur in the corresponding 
linear ODE's. 

4-1- Singularities and indicial exponents of the linear ODE of x^^^ 

We begin with the series for x^''''- for which we have obtained the linear ODE (modulo 
a prime) and have recognized all the singularities as given in (34). Our calculations 
provides a check on whether or not a diff-Pade analysis on a series too short to find 
the exact ODE can nevertheless yield enough information to locate the singularities 
precisely and determine the associated local exponents accurately. 

A diff-Pade analysis with just 400 coefficients (using approximating linear ODE's 
of order ten or eleven) already confirms, with 36 digit accuracy, the occurrence of the 
singularities given by the roots of the factors 

{l + w){l-3w + w^){l + 2w-4:w^) , (38) 
which are the Nickelian singularities labelled as 5), and the roots of the factors 

(1 - Tw + Sw^ -4w;3) (1 + 8^ + 20^2 + 15^3 +4^") , (39) 
with twelve digit accuracy. The roots of the factors 

(l-u;)(l +2w)(l-u; - 3^2+4 w^) (40) 

are confirmed with three or four digit accuracy, and, finally, the roots of the factor 
1 + Aw + 8 up arc obtained with just one digit accuracy. The roots of 1 + 3w + Aw"^ 
do not yet appear in the analysis using 400 coefficients. 

A generalized diff-Pade analysis is conducted by steadily increasing the order of 
the linear ODE and the degree of the polynomials, while looking for the roots which 
stabilise with increasing accuracy. For instance, using 1250 coefficients, the roots (39) 
now appear with fifteen digit accuracy, and similarly for the other roots, the accuracy 
increases. While the roots of the factor 1 + 3w + Aw'^ are not yet observed, we see 

^ Wc call difT-Padc the type of analysis detailed in [11]. For a given number of terms of the series, 
there is a linear ODE of order q that reproduces the first N — q terms but may fail for subsequent 
coefBcients. The same analysis, called "method of differential approximants" was described in [26]. 
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the emergence of a new polynomial not found among the set of singularities of : 
the singularity w = +1/2 is actually observed with twelve digit accuracy. 

Further increasing the degrees and the order of the hnear ODE to fully utilise 
the 1980 coefficients at our disposal, the accuracy is dramatically improved. The 
roots (39) are now confirmed up to 67 digits, the two complex roots of polynomial 
1 — w — 2>uP + Aw^ are observed with more than 51 correct digits, and the real root 
is observed with more than 7 correct digits. The roots of the factor 1 + 4m; + Sm;^ 
are now observed with 17 digit accuracy. The roots of (1 — w){l + 2w) arc seen with 
5 correct digits for w = 1 and 7 correct digits for w = —1/2. Finally, the roots of 
l + iw+Aw^ are seen with 4 digit accuracy, and the "new" (1 — 2iu) factor is observed 
with 27 correct digits. We summarise this discussion in Table 3. 



Table 3. This tabic shows the number of significant digits found, in the case of 
each singularity of x''''i from a differential approximant (diflF-Pade) analysis of 
series of length, respectively, 400, 1250 and 1980 terms. 



Label of singularity 


Associated polynomial 


400 


1250 


1980 




(l+w) 


36 






P{H) 


{l + 2w- 4ii;2) 


36 








36 






Pi^ 5i,4) 


{l + 8w + 20w;2 + 15w^ + 4w^) 


12 


15 


67 


^(^53.2) 


(1 - 7ii; + Bw"^ - 4w^) 


12 


15 


67 


^'('54,l) 
^^('53,2) 


(1 - w - 3w;2 + 4m;3) 


4 


15 


51 


(1 + 4m; + 8m;2) 




8 


17 


P(5 5/2 3) 
P(5 5/23) 


(1 + 2w) 


3 


5 


7 


{1-w) 


3 


3 


5 


J'('5/'3i,2) 


(1 + 3u; + 4«;2) 






4 


Unknown 


(1 - 2w) 




12 


27 



These diff-Pade calculations are, in fact, sufBciently robust to allow us to predict 
the minimum multiplicity of some singularities. When a given singularity (say w = 1) 
is observed, it is put in the head polynomial for a second run using more coefficients. 
If this singularity appears again it must be a double root (one exact root w = 1 and 
a second root w ~ 1 observed numerically with sufficient accuracy). Next the factor 
(1 — w)'^ is included in the head polynomial for another run and so on until no further 
occurrences of the given root are found. 

With the number of series coefficients at hand we find for the head polynomial of 
the ODE the following factors occurring with the indicated multiplicity: 

(1 +2wf ■ (1 -wf ■{! +Sw +4^2)2, 

while all other roots (except 1—16 w"^) occur with m.ultiplicity one. We are very close 
to the exact multiplicities of the exact linear ODE (sec (34)). 

This knowledge can then be used, in a kind of converging procedure, to improve 
the accuracy of our diff-Pade calculations. Having a totally unambiguous location of 
the singularities and a minimum value for the multiplicities, we revisit the diff-Pade 
calculations using the following Ansatz for the polynomial in front of the i*'' derivative 

w'-^ ■ (1 - 16«;2)'-2 . (^(1 + 2w) ■ (1 - u;))'""^* 

X (1 + 3w + 4«;2)'-''+2 . P(w)'-9+i . Q^(w) (41) 



Extreme series 



20 



with P{w) containing the other singularities with multiplicity one: 

P{w) = (1 + «;) (1 - 3«; + w^) {l + 2w- Aw^) (l + 4w + 8w'^) 
X (1 - 2w) {l-7w + 5u;^ - 4m;^) (l - «; - + 4w^) 

X (1 + 8u; + 20u;2 + 15^,3 + 4^;*) 

Here q is the order of the linear ODE and the index i = 0, ■ ■■ , q denotes the 
successive derivatives in the linear ODE (the actual exponents in the Ansatz are zero 
when a negative value is encountered). The Qi(w) are unknown polynomials, and 
Qq{w) at the highest derivative is included in order to handle the expected apparent 
singularities. 

The Ansatz (41) is then used in a diff-Pade analysis to obtain in floating point 
form, at each singularity, the associated critical exponents as roots of the indicial 
equation. Our findings are listed in Table 4. Let us explain how we display our 
results. Consider for instance the singularity 1 — w = 0. The successive roots of 
the indicial equation appear as the integers 0, 1, • • • , g — 5 (which is by construction) 
together with the integers 2, 3, 3, 4. The results for this case will be displayed as 
2,3^,4, but note that the roots 2 and 4 both appear twice, and the root 3 appears 
three times. Recall that the roots of the indicial equations appear in floating point 
form and we recognize these roots as being integers or half integers. The accuracy 
of the indicial exponents can be as low as three correct digits as is the case at the 
singularities given byl + 3w + 4«;^= 0. 



Table 4. Singularity exponent list for x'^^ , based on the exact {modulo a prime) 
linear ODE compared to those found by the (floating-point) diff-Pade analysis. 
The mod prime exponent list is followed by a number in bra<;es, which is the sum 
in (A. 7) of Appendix A. A final column gives the expected "normal" exponents 
using (E.l). 





Singularity 


Exponents (from 


Exponents 


Exponents 




Polynomial 


mod prime analysis) 


(diff-Pade) 


(E.l) 






l^ 2*, 3-*, 4^ 5^ 6^ 7^ 8^ 9^ 

10,12^ 15,25,(192) 










-2, -7/4, -3/2, -5/4, -l^ -1/2, 0*, 


-3/2,-1,0" 


-1 






1/2,1^,2^3,4,5,6,7, (549/2) 
-1,-1/2,0*, 1/2, l^ 3/2, 2^3^ 








{l + Aw^^ 


-l/2,o^l2 









4,5, (315/2) 
0^, 1*, 2^ 3^ 4^ 5^ 6, 7, 8, (56 + 91) 








l/t«i« 


o^l^2 






il + 2w)* 


2,5/2,3^,(41/2) 
2,3^4, (22) 


2,5/2,3^ 


3 




(l-wf 


2,3^4 


3 




{1 + 3w + 4w'^f 


0,1^2,(14) X 2 





1 




(1 + w) 


11,(12) 


11 


11 




(1 + 2w -4w^) 


11, (12) X 2 


11 


11 




(l-3w + w^) 


11, (12) X 2 


11 


11 


(l + 8u 




7, (8) 


7 


7 


(1 


-7w + 5w^ - 4w^) 
- w - 3w'^ + 4w^) 


5, (6) X 3 


5 


5 


(1 


7, (8) X 3 


7 


7 




{l + 4w + 8w^) 


5, (6) X 2 


5 


5 




(1 - 2w) 


7/2, (9/2) 


7/2 





Since we have obtained the exact linear ODE for x^^^ this case provides a valuable 
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test of our diff-Pade analysis. Furthermore, the Landau singularity analysis has been 
extended, in Appendix E, to include a power counting argument for the exponents at 
each singularity. All these results, the exponents from a diff-Pade analysis on 2000 
terms, the exponents from the (modulo a prime) linear ODE, and those derived in 
Appendix E are displayed together in Table 4. 

Even if we did not know the exact (modulo a prime) linear ODE of x^^\ our 
diff-Pade analysis can provide accurate information about the exact ODE. We note 
that the indicial exponents found above are accurate enough that we can be confident 
in their exact values, but the set of exponents may be incomplete. Let us consider the 
case 1 — 2w = to show what we mean. This factor was taken with a multiplicity 
of one in the head polynomial and the roots of the indicial equation show up as 
0, 1, • • • , g — 2, 7/2. Assume that in the exact linear ODE, the multiplicity is two. 
Then, the roots of the indicial equation will be 0, 1, • • • , q — 3, 7/2 plus an unknown 
exponent. If we had more series coefficients, with a further increase of the order and 
the degrees, this unknown indicial exponent may be obtained as another half integer, 
the exponent 7/2 or an integer. The dominant singular behavior at w = 1/2 would 
change accordingly. 

The comparison with the exact ODE results shows that the diff-Pade analysis 
on only 2000 terms (which are insufficient to encode the linear ODE for x'^^) is 
able to correctly give all the local exponents for the singularities together with the 
correct multiplicity. The singularity polynomial (1 -|- 3w -t- 4w'^) was used in the 
head polynomial of the linear ODE with a multiplicity of two instead of the correct 
multiplicity four. One should then obtain two local exponents. The missing local 
exponent was obtained as 0.87 instead of 1. 

These results give us considerable confidence that our numerical analysis of x^^^ in 
the following subsection, and of higher order susceptibility components subsequently, 
are completely correct. 

4-2. Singularities and indicial exponents of the linear ODE of x^^^ 

Similar calculations to those detailed in the previous section for x'^^ have been 
performed for x^^) in the variable x = w'^. 

In a diff-Pade analysis, increasing the order of the linear ODE and the degree 
of the polynomials, the singularities predicted by the ^^'^ model (labelled as 

6)P{'^ 64^2)P{'^ 65,1) and 4)) are obtained with increasing accuracy. 

We should note that the additional singularity 1 — 2w = 0, occurring for x'-^\ 
was seen in our diff-Pade analysis before we obtained the exact {modulo a prime) linear 
ODE. Let us detail for x'-^'' how our numerical procedure proceeds. First we check 
whether the "candidate" singularities appear as roots of the head polynomials, and 
we also check whether some of the other roots stabilize as the number of terms and 
the order of the linear ODE increases. If so this root is a true singularity. Next the 
well confirmed "candidate" singularities are put into the head polynomial and another 
run is carried out with more terms to confirm the "new singularity" . 

Table 5 shows the results of three specific runs. The third column shows results 
from a run using 387 terms with order 12 ODE. An unknown (with respect to our 
candidates) singularity is seen to stabilize around a value consistent with 1 — 8x = 
to an accuracy of 3 digits. Keeping the same number of terms but increasing the order 
to 16 this root would change if it were a root of the apparent polynomial. The results 
in the fourth column show that the new unknown singularity remains and is seen with 
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4 digit accuracy with respect to its exact value. Finally all the other singularities are 
put into the head polynomial and the fifth column shows the results for an order 31 
analysis using 997 terms. These calculations show the existence of a new singularity 
at X = 1/8 and thus arising from the factor 1 — Sw"^ = in the head polynomial of 
the true ODE. 

Table 5. This table shows the number of significant digits found, in the case of 
each singularity of x^^\ from a diflf-Pade analysis (a; = w"^). 



Singularity 


Singularity 


387 terms, 


387 terms, 


997 terms. 


Label 


Polynomial 


order 12 


order 16 


order 31 


P('6) 

P{^6/H) 

P(26) 


1-x 

1 - Ax 


28 

30 


29 

30 




l-9x 


30 


30 




^(^64,2) 


1 - 25a; 


13 


14 




^(^64,2) 


1-X + 16a;2 


8 


10 






1-10X + 2Qx^ 


10 


12 




Unknown 


l-d>x 


3 


4 


26 



Various runs (with linear ODE's ranging in order from 26 to 31) show that the 
factor (1 — 4 a;) occurs with a power five which leads us to the following Ansatz: 

a;'-i • (1 - 16 x)'-^ • (1 - 4 a;)*-«+^ • P(a;)'-«+^ ■ Qi{x) (42) 

for the polynomials in front of the i*'' derivative in the linear ODE where P{x) reads: 

P{x) = (l-a;)(l-9x)(l-25x)(l-a;-16x2) 

X (1-10.t + 29.x2)(1-8.t). 

The singularities and corresponding exponents for x'^-' are summarised in Table 6, 
and we also show the exponents predicted by the local analysis given in Appendix E. 

Table 6. Singularity exponent list for x'^^ found by the (floating-point) diff-Pade 

analysis together with the expected "normal" exponents (E.l). 



Singularity 


Singularity 


Exponents 


Exponents 


Label 


Polynomial 


(diff-Pade) 


(E.l) 




X 


0,-1^-1/2 






1 - 16a; 


-3/2,-1,05,1 
-12,0^,-1/22,1/26 


-1 




1/a; 




P(26) 

P(-"^6/24) 

P(26) 


\-x 


33/2 

11/2,13/22,15/2,33/2 


33/2 


I -Ax 


13/2 


l-Qx 


33/2 


33/2 


^('64,2) 


1 - 25.T 


17/2 


17/2 


^('64,2) 
^('65,1) 


1 — a; + 16a;^ 


17/2 


17/2 


1 - lOa; + 29a;2 


23/2 


23/2 


Unknown 


l-8a; 


7 





Again we see that with some 1600 terms (which are insuflacient to encode the linear 
ODE for x^6^) the diff-Pade analysis is able to give the correct exponents in agreement 
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with Appendix E. Depending on the multiphcity of the singularity polynomial put into 

the linear ODE, the exponents may be incomplete. Again we find rational exponents 
for all the singularities of the linear ODE for x'^^K 

From our diff-Pade calculations we thus have confirmation that the singularities 
encoded in the head polynomial of tlic^ as yet unknown linear ODE for x^^^ are (at 
least) the ones of together with the roots of a new polynomial 1 — 8 w^. Note that 
these 1 — = additional singularities lie on the unit circle \s\ = 1. At this stage, 
and since all the series coefficients of x^^-* have been used, one may ask if there are 
other singularities like w"^ = 1/8 "still to be discovered"? Without the exact linear 
ODE of x^^^ we cannot give a definitive answer to this question. 

5. Diff-Pade analysis for x^"'\ n > 7 

It is known [8] that the first non-zero coefficients in x^"^ appear at order in w. 
The high- or low-temperature series for x up to iV = 2000 then contains contributions 
from all the odd, respectively even, x*-"^ up to n = ^/N. 

Since we know the first 2000 series coefficients for x*'^'' and the first 1630 
coefficients for x'^'' a-s well as the series for x'"\ n < 4, up to an arbitrary number 
of coefficients, one may ask whether the total x with these lower x*-"' terms removed, 
can yield any information about the singularities that should occur in the linear ODE 
of x^"\ n > 7? As far as x^^'' and x^^^ are concerned and in view of the "limited" 
analysis done here for x^") , n>7, our conclusion is that there are no new singularities 
that are not in the "known" set (i.e. irrespective of the index n). 

5.1. High temperature analysis 

The diff-Pade analysis of the long high-temperature series 

XH - X^^' - X*^^ - X*'"^^ = 128 w^^ + 25088 w'^^ +■■■ (43) 

amounts to looking for the singularities of a "linear ODE" approximation to the 
(infinite non-holonomic) sum, x'^'' +x'^'' +x'^^^'' +• • • , as roots of the head polynomial 
of a given diff-Pade approximant. These singularities occur grouped together (that 
is to say we find several singularities simultaneously) and each singularity can be 
attributed to a given x^") according to the polynomials given in Appendix C, i.e. 
arising from the (f^-* model. 

With an order eight linear ODE and 312 terms, the roots of the following 
polynomials (besides 1 — 16w^ = 0) are recognized: 

Pi : 1 + 2u; - 8u;^ - 8w^ = 0, 

P2 : 1 + 2w - w'^ - = 0. 
P3 : 1 - 5w + 6w^ - = 0. 

These polynomials can be identified with all the Nickelian singularities corresponding 
to x^^). The roots of the polynomial Pi are obtained with 5, 10 and 12 correct digits, 
while two roots of the polynomial P2 are obtained with 3 and 9 correct digits. 

A summary of various runs using 500, 900 and 1956 terms is given in Table 
7, where we display the number of significant digits found, in the case of each 
singularity of the partial high-temperature susceptibility (43). The Ansatz for the 
linear ODE search is to include only the factor (1 — 16w^)*, where i is a positive 
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Table 7. This table shows the number of significant digits found for each 
singularity of the partial high-temperature susceptibility (43) using a diff-Pade 
analysis. The singularities equal those predicted from 3"^' as occurring in x^"^- 
The second column gives the polynomial factor, the zeros of which corresponds to 
the singularity locations. Subsequent pairs of columns give, firstly, the number of 
significant digits found for that singularity, and secondly the order of the linear 
ODE for which this singularity was found. 







500 terms 


900 terms 


1956 terms 


n 


polynomial 


Digits 


Order 


Digits 


Order 


Digits 


Order 


7 


1 — 5w + — 




18 


12 


33 


19 


50 


18 


7 


1 + 2w - u? - 




16 


11 


28 


15 


58 


18 


7 


1 + 2ir - Sir- - Sir' 




lo 


1 1 

1 i 




1 

1-) 


■)S 


lo 


y 


1 — w 












15 


14 


9 


l + 2w 








6 


13 


30 


16 


9 


l + 3w- up 




5 


11 


16 


15 


45 


15 


9 


1 — 6w + 9w^ — tu^ 




7 


12 


20 


13 


45 


18 


9 


1 - ;W - 








8 


15 


30 


18 


9 


1 - 12(c- + 8«''' 








14 


14 


40 


20 


11 


1 - 9w + 2Mi^ - 35w^ + 15^"^ - 


5 

w 






7 


13 


30 


18 


11 


1 + 2w - - 2w^ + Aw"^ - It 


5 










23 


17 


11 


1 + 2ir - !()/(•- - 24(r'"' + -iSir'- + ? 


2 ((■""' 










20 


20 


13 


1 - llw;H h 'u/* 












16 


16 


13 


l + 2w- 20w;2 _^ g4^6 












7 


18 


15 


1 + 2w - Aw^ 












5 


14 


15 


1 - 9«; H h 












5 


14 


7 


1 + 12m; + Mw^ H h Aw'^ 








6 


13 


12 


17 


7 


1 - Su; - lOu;^ H 16w^ 












4 


12 


9 


1 + 16«; + 104«;^ H h 4u;** 










7 


17 



integer discussed in the previous section. The first block of singularities correspond 
to NickcUan singularities, given as Case 2 in Appendix C. 

Again we have confirmation that some of the singularities of the linear ODE 
of are actually singularities of the linear ODE of the x^"''t- Prom Table 7 we 
note that the first singularity polynomials to appear are the Nickelian ones. All the 
Nickelian singularities for n = 7 and n = 9 are confirmed. Of the four singularity 
polynomials for n = 11, one is missing and among the five singularity polynomials for 
n = 13, three are missing. We also note that once the Nickelian singularities for a 
given n have appeared the non-Nickelian ones begin to show up. 

t Note however that without an "exact" series for each x^") , one cannot safely attribute, for instance, 
the non-Nickelian singularities polynomial 

(1 + 16w + 104«)2 + 352^3 + mOw* + 672w^ + 336w^ + 63w'^ + 4w^) 
to x'^^ ■ These singularities may well come from x^'^^ , but not be predicted by . 
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5.2. Low temperature analysis 

Similarly, the diff-Pade analysis of the long low-temperature series {x = w'^) 

XL - X^^^ - X^^^ - X^^^ = 256a;^2 + 65536a;^^ + • • • (44) 

amounts to looking for the singularities of the "linear ODE" of the (infinite) sum, 
-j^(8) _|_T^(io) _|_T^(i2) _|_. . . j-oots of the head polynomial of the diff-Padc approximant. 

The Ansatz for the ODE search is to include only the factor (1 — 16a;)*. Here we 
show the results in Table 8 for two orders fully utilising the series coefficients at our 
disposal. The first block of singularities correspond to Nickelian singularities given as 
Case 2 in Appendix C. 

Table 8. The number of significant digits found for each singularity of tfic partial 
low-temperature susceptibility (44) using a diff-Pade analysis. The singularities 
are those predicted from 3"^' to occur in x'"). The second column gives the 
polynomial factor whose roots give the singularity locations. Subsequent columns 
give the number of significant digits found for that singularity at orders 14 and 
20. 



n 


polynomial 


Order 14 


Order 20 


8, 12, 16 






l-4x 


69 


72 


8, 16 






1 - 2a; 


51 


52 


8, 16 






1 - 8a; 


74 


75 


8, 16 






1 - 12a; + 4a;2 


73 


79 


10, 12 






1-x 


21 


24 


10 






1 - 5.T 


43 


46 


10 






l-7x + x^ 


42 


44 


10 






1 - 12a; -1- 16a;2 


49 


53 


10 






1 - 15a; + 25a;2 


55 


59 


12 






1 - 9x 


12 


16 


12 






l-3x 


10 


12 


12 






1 - I2x 


9 


13 


12 
12 






1 - 14,T + .T^ 

1 - 8.r + 4,i:^ 


35 
(i 


40 
() 


14 




1 


- 21a; + 98a;^ - 49a;^ 


21 


24 


16 


1 - 


- 24a; + 148a;^ - 176a;^ + 4a;4 


10 


12 


8 




1 


- 20x -1- 16a;^ - 16a;=* 


8 


7 


8 


1 - 26a; ^ 


- 242a;2 - 960x^ + 1685a;'' - llSSa;^ 


6 


6 


10 




1 - 


- 24a; + 128x^ - 289a;^ 




4 


10 


1 - 


46a; -|- 866a;2 -|- 56642a;9 




6 



As for the high temperature analysis we have confirmation that some of the 

(n) 

singularities of the linear ODE of arc actually singularities of the linear ODE of 
the x*"""*- All tli6 Nickelian singularities for n = 8, 10, 12 are confirmed. Of the five 
(nine) Nickelian singularity polynomials for n = 14 (n = 16) four (eight) are missing. 

As for the new singularities, 1 — 2w — for the ODE of x'-^^^ and 1 — 8uP = for 
the ODE of x^^-*, we should say that these singularities occur for higher index n. Thus, 
as far as x^^^ and x^^^ are concerned, and in view of the "limited" analysis carried 
out for x^"\ n > 6, there is no new singularity discovered that is not in the "known" 
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set (given by the Nickelian singularities and the singularities of the ^f^' integrals, 
irrespective of the index n > 6). 

5.3. Local exponents 

We turn now to the indicial exponents at some of the singularities found in the analysis 
of (43) and (44). 

Here, the situation may seem different from the equivalent analysis made for x*-^-* 
and x^^^ . For instance, when we consider the indicial exponents for the singularities 
that should occur in the still unknown linear ODE for, e.g. x^^\ the series that we 
analyse contain contributions from all the even x'"'' up to VTV. For instance, for 
(1 — 4x) = 0, this singularity is also a singularity of the linear ODE of x^^^\ x'^^\ 
etc. One may thus expect to obtain for the Nickelian singularity (1 — 4x) = the local 
exponent [4] (n^ — 3)/2, for n = 8, 10 and n = 12. The differences between all three 
exponents being integers, it is the value 61/2 corresponding to n = 8 which should 
appear. What we have obtained is indeed that for the Nickelian singularities, when 
the exponents have stabilized, they agree with (n^ — 3)/2, where n is the lowest index. 

Here we give some other examples of local exponents. For the non-Nickelian 
singularity polynomial (1 — 26x + 242a;^ — 960x^ + 1685x^ — 1138a;^) corresponding to 
n = 8, the local exponent is 23.5 (with 4 digits for g = 14, 5 digits for q = 17 and 3 
digits for q = 20). For the non-Nickelian singularity polynomial (1 — 20a:-|-16a::^ — 16a;^) 
corresponding to n = 8, the local exponent is 18.5 (with 2 digits for q = 14, 3 digits 
for g = 17 and 2 digits for q = 20). Both exponents agree with (E.l) with n = 8 and 
respectively m = 1 and m = 2. 

Similarly for the high temperature analysis, and for instance, for the Nickelian 
singularities given by roots of the polynomials Pi , P2 and P3 corresponding to x^^^ , 
the indicial exponent 23 appears with 10 correct digits using 480 terms at order 25. 
The exponents for the non-Nickelian singularity polynomial (l-|-12w-|-54'u;^-|-112w^-|- 
105w^ + 35w^ + 4:W^), which should correspond to x^^\ shows up as 17 with 2 correct 
digits. 

Let us close this analysis with the following remark: the diff-Pade analysis of the 
sum (43) has shown the singularity w = —1/2 which was attributed to x'^^''- Actually 
w = —1/2 is known to be a Nickelian singularity for x^^^- Prom the analysis of the 
we found that w = —1/2 is also a singularity of $2^- So this singularity may 

well be attributable to x^^-* • The diff-Pade analysis gives for this singularity 39 as the 
local exponent, which is the value predicted by (n^ — 3)/2 for x'^''§- 

5.4. "Indicial" exponents of non-holonomic sums 

We showed in previous subsections that even if the series of x^^) is "polluted" by the 
terms of the other x''"-* the diff-Pade analysis is efficient enough to give the correct 
singularities and the corresponding indicial exponents. 

Recall that the sums (43, 44) arc believed to be non-holonomic, thus our linear 
ODE search (diff-Pade analysis) is just an approximation and another way to encode, 
via the singularities and their local exponents, the information contained in the series 
coefficients. 

§ We may thus conclude that either the Unear ODE of x^^' docs not have this singularity or it has 
this singularity, but the exponent is compatible with 39, i.e. it may be {p)P with p > 39, p > 1. 
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For this series with infinitely many singularities, generically not suitable for a 

difF-Pade analysis, it is interesting to see how the local exponents actually appear for 
those singularities occurring in all x''"'') i-^- w = ±1/4. 

To more fully appreciate the results of this section, consider the "holonomic" sum 
■^(1) +;^('^) +X^^^ ■ One may ask what arc the local exponents that appear for w = 1/4, 
knowing the exponents for each term? Recall that at w = 1/4, the local exponent for 
the linear ODE of x*^^ is -1, for x^^^ they are -3/2, (-1)^ and (0)^ and for the linear 
ODE of x'^''^ (from diff-Padc) the observed exponents are —3/2, —1 and (0)"*. 

A diff-Pade analysis on the holonomic sum x'-^'' + X^'^^ + X^^^ should show all 
these exponents (with possibly a change in the multiplicity due to auto-cancellation) 
and possibly other exponents that differ by integer values from the exponents in each 
individual x'-"'' term. This last possibility comes from a cancellation between the 
initial terms of the series. 

Putting the known singularities into the head polynomial of an order fifteen linear 
ODE and using 860 terms we obtain: 

w = 1/4, -3/2, -1, 0^, 1. (45) 

Thus we see that here is no new exponent not differing by integer values from the 
exponents of each individual x^"^ term. 

Let us now return to considering the non-holonomic infinite sum, xh = x^'^'^ + 
-j^(3) + . . for which a diff-Pade analysis (with some 2000 terms and orders 15 and 
18) gives, at w = 1/4, the local exponents: 

w = 1/4, -1, -1/8, (3/8)2, (i5/g)2^ (35/8)2 (^g) 

Here wc sec the generation of new exponents that can only be seen as a consequence 
of the non-holonomic character of the sum. Note that the new exponents are still 
rational numbers. Note also that the new exponents appear not only in a diff-Pade 
analysis of xh, but also in xh — x'^^^ — X^''^ or in Xh — x'^^ ^ X^''^ ^ X*^'''''- One has 
some kind of "self-similarity" in x as far as the indicial exponents are concerned. Note 
that the multiplicities above may increase with more terms. 

The local exponents are those of the full xh around w = 1/4. They can be 
checked in 

Cl _ s''')i/^ 

— Xh = const ■t-'^'^F+ + Bf (47) 

given by Orrick et al, using and Bf in the Appendix in [8], by switching from the 
variable r to the variable {w — 1/4). 

The following section presents a simple model showing the mechanism by which 

a resummation of the infinite number of logarithmic singularities proliferating in the 
holonomic n-fold integrals x*-"'' can lead to the known power-law singularities in x- 

6. Convergence of the n-particle sequence 

It was already observed by Wu et al. [2] that for the leading divergence of the 
susceptibility proportional to |t|^'''/'* (where r = (1/s — s)/2), the sequence of 
partial sums of either x^^"^ or ;^(2n+i) appears to converge exceptionally rapidly to 
the corresponding susceptibility below or above the critical temperature. This has 
been confirmed numerically to higher order by Bailey et al. [25], who estimated that, 
asymptotically, the amplitudes /„ of the leading divergence of x^"^ are in the ratio 
In+2/In < 1/1000. On the other hand, the leading correction terms in x^^^ and x^^) 



Extreme series 



28 



are easily seen to be of order |t|^/^ and |rp/'' log |r| respectively whereas the leading 
correction in x approaches a constant in the limit r — > 0. In fact the exact solutions 
for x'^' and x^^\ found by Zenine et al. [10, 12] make it plausible that every partial 
sum of x^^"^ or x^^""*"^^ will, asymptotically in the limit r — > 0, give a vanishingly 
small contribution to the leading correction term. This raises the issue of the nature 
of the convergence of the n-particle sequence to the susceptibility. 

With the exceptionally long series now available for x^^^ and x'^^ we can, even 
in the absence of an exact linear ODE solution, make numerically precise estimates 
of the correction terms in these functions. A scheme that works well is a combination 
of unwanted singularity suppression as described in [3] and function fitting. The 
latter is done by assuming various combinations of powers and powers of logarithms 
with unknown coefficients and generating the corresponding series in s. The unknown 
coefficients can then be estimated by least squares fitting to the highest order 
terms of the exact (unwanted singularity suppressed) scries. This fitting is typically 
done iteratively with the leading order exact coefficient values substituted for the 
numerically estimated ones from an earlier iteration. 

Then, by combining this information with what is already known for the lower 
orders x^"^ > we can build a plausible "toy model" Xtoj/ ^^^^ '^^^ ^® easily extrapolated 
to n — > 00. This gives us, at least qualitatively and even semi-quantitatively, a 
picture of the convergence of the partial sums of x^"^ to the full X- The details of 
this extrapolation procedure via Xroy described below. 

6.1. Behaviour of x^^^ 

We begin by providing a summary of the behaviour of x''"'*'"^'' near the ferromagnetic 
critical point. 

The normalization factor included for convenience on the left-hand side of the 
equations below is Air^nl while the factor i/i is crucial to give series increasing in even 
powers of t only. An empirical observation is that t = (1/s — s)/2 appears almost 
universally in all formulae in the combination |t|/4 so we have defined T4 = |t|/4. The 
known x^"-* terms arc| 

^VS-X'"' = 4- Qr,-^/^ +rr) -r^ + 0(rn, 

4-^/i•x(^) = ^•Qr-/*+rr) + (logr4 + §)-rl/nO(rr), 

- (log^ + f + y + i) • -1'' + 

+ (2iT^ + ^) logr4 + 34.3411462895 . . • tI'^ + 0{tI'^-% 
96n'V-s-X^'^ = 967+7r^. Qr-^/^ + r]/^) - (^log V4 + | log V4 

II The formula given here for x'*^ corrects misprints in [15]. 
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671-2 _^ . ^Qg2 _^ 305.6550541085 . . . logT4 
+ 375.271992213596336733341995793^ t^^ + 0{t^^^~^), 
mn'VS-X^'^ = 480/e-7r5. Qr^/^ + r]/^) + (^log^ ^4 + ^ logV4 

40n ^ 2275 \ ^^^3 _^ 1437.2558956691 . . .In^ T4 
3 1° / 
+4238.8509988858869410798745 logT4 

+ 5284.244417602341195112209^ • t^^* + 0{t^^*^^) (48) 
while the complete summation of all x^") yields the exact 

^/s■x= 1.00081526044021... ^7-4"^^^ +^1^*^ - 0.1041332450938... 

+(0.0323522684773 . . . log(r) - 0.074368869753 . . .) • r + 0{t^) (49) 
for r > and 

^ 1.00096032872526... /I .7/4 imA „ ,„.,oo^..^^o„ 
V^-X= TT^ 0^4 +^4 - 0.1041332450938... 

+ (0.0323522684773 . . . log(-T) - 0.074368869753 . . .) • t + 0{t^) (50) 

for r < 0. We have used 0(t4^^~'^) to indicate 0{tI^^) with logarithmic corrections. 
We also distinguish between numerical constants. Those with trailing . . . are known 
to much higher accuracy but have not yet been recognized in terms of elementary 
constantsf. The others are believed, but not guaranteed, to be accurate to the number 
of digits given. 

An important remark about the form of equations (48) is that we have split the 
tI^^ contribution to create the combination (^r^ '^^^ + t^^^) which is the scahng part 
of the full susceptibihty in (50). The original motivation for this was the observation 
by Orrick et al. [8] that the terms in the expansion of the scaling function in the 
susceptibility are '^/4+2" without logarithms, while in the expansion of the background 
part of the susceptibility there are only integer powers of T4 but now with logarithmic 
corrections. A simple way to reproduce the scaling part of the susceptibility from the 
X^"^ sum is to replace every r^"^^^ by the scaling function, as we have done. Then the 
sum rule on the J„ that ensures the correct r^^^'' amplitude in the susceptibility will 
also automatically yield the correct scaling function. Such separation into "scaling" 
and "background" in individual x^") would appear to be completely arbitrary except 
that we now find that the "background" remainder in x^"^ has a simple dependence on 
n in which the same formulae apply simultaneously to both odd and even n. 

This was quite unexpected and the separation may well only apply to the leading 
correction term we are concerned with here. As a practical matter, it leads us directly 

t The exact values for /g" and have been given by Tracy [27] (see also note added in 
proof in [15]). Highly accurate numerical values for and 1^7 are given by Bailey et al. [25]. 
For the constants in x sec Orrick et al. [8]. The constant of x^^' in (48), to higher accuracy is 
34.3411462895318287823760337015134958201194891883303683732355064700755161069888189217774 
8221462488591136755297844305262613713337765 • • • 
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to the conjecture 



-(-1)" 



log" T4+n{n-^]- log" ^ T4 



/n\ / 2 (27r2-19)n 125 - 127r2\ , „_2 

/nX /3^ (4.--37K ^ (374-37.> _ 0.4414526IO. . .) • log"- 



1/4 



+ 0(r9/*-^) (51) 

based on the known low order results and applicable equally to the logarithmic and 
non- logarithmic terms. Note that the numerical constant in (51) is simply related 
to the constant in x^^^ deduced from its linear ODE solution. Specifically, one has 
-0.44145 . . . = 34.341 . . . + tt'^/A - 149/4. 



6. 2. Resummation of the toy m,odel 



The conjectured general term (51) forms the basis for our x'xay^^- ^he leading 
divergence in x'"^^') now combined with the scaling function, is not under 
consideration here since wc have nothing to add to what is already known. Of the 
correction terms in (51) we will capture exactly the leading log" T4 into our toy model. 
The appearance of binomial coefficients in the next three terms is suggestive of a 
formula like (log T4 + n)" which would capture the leading n dependence correctly. 
We do not have enough information to be confident about the behaviour of any lower 
order logarithmic terms and so at this point our model becomes dictated by the criteria 
of simplicity. The numerical values of the lower order terms suggest the simple formula 
(logr4 + An)" as reasonable, where A « 1 is some as yet undetermined constant. If A 
is very close to 1.0 then our constant term (An)" is probably an underestimate but it 
is worth remarking that in the t ^ limit, logr4 is negative. This implies successive 
terms in the expansion of (logT4 + An)" alternate in sign and the final value, just as 
in the exact x("+^\ is the result of large cancellations between terms. Thus it is more 
important that we capture correctly the "smoothness" by which successive terms vary 
and this is hard to estimate from the limited data available. Let these caveats be 
understood. Let us denote by A the correction to scaling. Then our toy model for 
the correction to scaling terms is 

.1/4 

4 

^ 47r"n! 

V^-AxToy - ^^/i•Ax^'4 (52) 

with the sum understood to be over even n for T < and odd n for T > Tc- 
As T — > 0, the sum in (52) is dominated by large n and can be replaced by an 
integral that is easily treated by steepest descent methods. The value of the integral 
will depend on the undetermined A and we will choose A such that, in the limit 
T — > 0, "v/s • Axtoj/ ^Toy is a constant. That this is possible is verified by explicit 



Vs-Ax'^;'^ ^ --^•(log(l/r4) -An)^ 
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calculation below. The reasonableness of our toy model can then be judged by how 
close Axoy is to the exact A — —0.1041 . . . from (50). 

The dominant n dependence of the integrand in the integral approximation for 
^/s ■ AxToy is the factor exp[nlog(log(l/T4) — An) — n log(n7r) + n] and this has a 
maximum at rip where: 

np= • log(l/T-4), p • A + log(7rp) = 0. (53) 

The exponential at its maximum evaluates to exp(plog(l/r4)) = l/rf which 
establishes the dominant r dependence of the integral as a function of A since p = p{X) 
is the solution of the transcendental equation in (53). Expanding the exponential 
about its maximum then gives the approximate ^/s ■ Ax^^^^ and the required sum 
in (52), namely 

^•^^-^ " ^yrfnexp[-(l+,A)^(n-MV(2M] = (54) 

valid in the limit r ^ both above and below Tq. The choice p = 1/4 is now seen 
as necessary and we get from (53) and (54) 

A = 4 log(4/7r) « 0.966, 

Aroy = lim AxTo. = -J^^^^^^^ - "O-IOOT. (55) 

The close agreement of the asymptotic amplitude with the exact A k, —0.1041 gives us 
confidence that (52), with the specific choice of X from (55), will be usefully predictive 
for finite t and n. 

The result of numerical computation for finite r and n is shown for T < Tc in 
Fig. 1. A similar plot could be made for T > Tc- The agreement between exact and 
toy partial sums at low order is of course by design and the agreement in the limit 
T — > has already been remarked on in connection with (55). The striking feature of 
the shifts in log^oC^''') with order can be deduced from (53). The Up in (53) which 
corresponds to the n of maximal contribution is also the n characterising the transition 
region of the partial sum approximations to ^/s ■ AxToy between the asymptotes 
and Axoy = —0.1007. Then, since each additional order is a change in n of 2, we 
deduce a transition region shift of Alog(l/r4) = 8 • (1 + log(4/7r)) « 10, or about 
four decades in temperature. 



7. Natural boundary: power spectrum analysis 

In this section we are using the power spectrum analysis introduced in Sec 6.3 of 
Orrick et al. [8], though not fully described there. The model of n-fold integrals [14] 

(n) 

has shown that the singularities of the linear ODE of occur in the linear ODE 

of the higher order . These model integrals are assumed to "mimic" correctly 

the x^"' in terms of the locus of the singularities. With the singularities obtained 
here from diff-Padc analysis of x^^^ and x^^^ we have confirmation of this feature. 
If we focus on the singularities which are on the unit circle \s\ = 1, the Nickelian 
singularities of x^^^ > (i-e. {l — w){l + 2w) = 0) are also singularities of the linear ODE 
of x^5). 

However, the detailed Landau conditions analysis (see Appendix B) proves that 
the above singularities lying on |s| = 1 do not occur on the principal disc of X^^^- 
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Figure 1. The contribution to the "background" part of the susceptibiUty from 
the partial sums 52 = x'^^^Vs, S4 = S2 + x'-'^^Vs, 86=84 + x'^'v^ and the 
exact xVs based on the formulae (1) and (2) exclusive of the scaling ^^''+r^^'* 
combination. Also shown are the corresponding toy model Sn based on (52) now 
extending to larger n. 

Indeed, this is just a special case of the general theorem that applies to all the Case 3-5 
singularities. Furthermore, the analysis in Appendix D suggests that the singularities 
associated with the ODE head polynomial factors (1 — 2w) for x^'"'^ and (1 — 8w^) 
for x^^) are not even Landau singularities of the integrals x^^' a^^d x^^\ at least in 
a Landau singularities analysis that does not fully take into account the Fermionic 
term. 

Here we confirm the limited result of the absence of these singularities on |s| = 1 
on the principal disc using a technique [8] based on the fast-Fourier transform of the 
series ofx^"^- The scries is first appropriately modified by removing the dominant 
singularities and smoothing. 

The Nickelian singularities sq have very small amplitude, roughly (1 — s/sq)" 
for x^"-*. For x*-^-* the effect on the scries coefficients (for series in the s variable) 
at iV = 2000 is roughly 1/(2000)^3 = lO'^^^ times smaller than that from the 
ferromagnetic singularity. The amplitude of the dominant ferromagnetic divergence 
being known this contribution can be subtracted. As the smoothing continues (see 
details in [3]), the series coefficients will start to decay rapidly with order N and this 
must be corrected for by multiplying, at each step, by some power of A'' to again make 
the series coefficients roughly constant in A^. Once the series have been smoothed and 
multiplied up by about A^^^ the FFT is performed. 

The FFT is obtained by using say 512 coefficients a„+jv starting from some large 
N to utilise the highest order coefficients available. We get 

bm = ^a„+jvexpr2 7rin^j , m = 1, 2, ... ,512. (56) 
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Figure 2. Power spectrum of x^^' series on logj^Q scale vs "frequency" m 

and the power spectrum = bm ■ &5i2-m which runs over 256 points. The 
"frequency", m, is directly interprctablc as an angle 9 = 7rm/256 giving the 
locations of the singularities exp(±i0) in the complex s plane. Windowing a„ 
an sin^(n7r/512) can be used to reduce background (this broadens the spectral lines). 

The power spectrum of x'^^'^ series is shown on Fig. 2. The curves correspond to 
various values of the integers p = 2, • • • , 5. From left to right, the first and third 
spikes correspond to the singularities of 1 — 3w + w'^ = 0. The second and fifth spikes 
correspond to 1 + 2w — = and the singularity w = —1 appears as the fourth 
spike. 

The singularities w = 1/2, w = 1 and w = —1/2 are not seen. Their positions 
are indicated, respectively, by the vertical lines at the top of the spectrum. 

Similar analysis of the series of x^^^ has shown the non-occurrence of the spike 
corresponding to the non-Nickelian singularity = 1/8 which is on |s| = 1. 

The FFT analysis does show, in all examples so far, that there is no evidence 
for any singularities on the |s| = 1 boundary of the physical sheet other than the 
Nickelian singularities. Thus the possibility that there will be "destruction" of the 
natural boundary by cancellation of the singularities accumulating on the unit circle is 
becoming more and more remote. 

8. Conclusion 

We have calculated series of some 2000 or more terms for the magnetic susceptibility 

(high- and low-temperature regime of the isotropic Ising model) as well as for the 
individual five- and six-particle contributions x'"^^ and x*^^-' . This was achieved by use 
of modular arithmetic which amounts to generating the coefficients modulo a set of 
primes, then using the Chinese remainder to obtain the exact coefficients. 

By extending the series for x^^^ modulo a single prime to 10000 terms we have 
discovered the exact linear ODE for x^^' modulo that prime. The Russian-doll 
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structure previously found to occur for x^^^ a^^d x^^\ and conjectured for the linear 
differential operators of the x^"^ , is actually confirmed for x^^^ ) weZZ as a stronger 
direct sum structure. 

We have described our present algorithm for finding the linear ODE satisfied by 
a scries S{x) in some detail because it is robust and particularly efficient. It is based 
on an ansatz of a rectangular array of coefficients of size fixed by degree D in x and 
order M in a;d/da;. The use of the operator xd/dx rather than d/dx guarantees that 
the resulting ODE is Fuchsian. Now suppose we have found the ODE with minimum 
D for given M. If is the number of series terms required for this D, M combination 
then an empirical observation is that there exists a linear relationship between N, M 
and D that enables us to infer the order Mq of the minimum order ODE required 
by >S'(a;) without having to obtain the minimum order ODE itself or even having to 
obtain any singularity information whatsoever. We have no explanation for this linear 
relationship and leave it as a challenge for the reader to provide an analytic basis for 
our "magic formula" (30). 

The exact long series have been used in diff-Pade analysis to obtain the 
singularities that should occur in the linear ODE of x^^^ and x^^^. We have 
confirmation that the singularities of the ODEs of the n-fold integrals and 

^^fp (which differ from x*-'"''' and x''^'' by the absence of the Fermionic factor) are 
all singularities of the ODEs of x*-^-* and x^^-* • However, our diff-Pade analysis showed 
that these linear ODEs have additional singularities, namely w = 1/2 for the linear 
ODE of x^^^ and vP' = 1/8 for the ODE of x^^'; ''^ot predicted by the corresponding 
We see that the Fermionic factor does not affect the singularities corresponding 
to but it does bring additional singularities to the linear ODE. From a Landau 

singularity analysis viewpoint this is not a surprise: A new factor cannot move a 
singularity, but it can either give rise to cancellations or introduce new singularities, 
and we have seen that cancellation does not occur. Our exact mod prime ODE for 
X^^^ confirms the accuracy of the diff-Pade results and in particular proves that the 
ODE for x^^^ carries the extra singularity at w = 1/2. 

We do not know whether or not these extra singularities of the x^"'' ODEs are 
singularities of the integrals (5). The most common case is that in which all the 
singularities of the ODE and the integral are the same; this is the situation for x^^' 
if one takes into account its various analytical continuations (cf. first footnote in 
Appendix B). As a specific example, both the Landau analysis and the ODE for x^^^ 
predict a singularity at s = (— l + i\/7)/4 for which it is to be noted |s| < 1. While x^^^ 
on the principal disc is not singular at this point there exists an analytic continuation 
of x^^) that is|. Our analysis of toy integrals in Appendix D provides an example of a 
difference between the singularities of the ODE and those of the integral. The ODE 
for the toy analog of x^^'' has singularities at w; = 1/2 and uP' = 1/8 but our Landau 
analysis of the toy integral fails to find singularities at these points. This may be a 
genuine distinction, or it may be that we missed something in the Landau analysis, 
or that the Landau analysis can't bo guaranteed to give all singularities. It would 
be extremely useful to have other (simpler) examples where it can be shown that the 
ODE defined by an integral has more singularities than the integral. 

A further check on the occurrence of the singularities (of the linear ODE's of 
$^^) has been made for some x^"\ i^ ^ 7, by subtracting, from the full Xi the 

X This can be inferred from the connection matrices provided in [15]. 
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known x^^^ x^^^ and the now long x^^^ series. The diff-Pade analysis again confirms 

(n) 

a large number of singularities occurring in the linear ODE's of the that can, 
thus, be attributed to the corresponding linear ODE of x^"^. A similar analysis has 
been performed for the low- temperature regime. 

While those diff-Padc analyses only yield approximate linear ODE's, they are 
efficient enough to give the indicial exponents at each singularity with good accuracy 
in most cases. An example of failure is that in the ODE for x^^^ at w = 1/4 where 
we detected —3/2 but not —2,-7/4 and —5/4, exponents also known to be present 
from our exact mod prime analysis. Thus while in all cases we have observed that 
indicial exponents are rational numbers, in the absence of this example we might have 
erroneously concluded that indicial exponents are either integer or half-integer. As a 
particularly striking success we obtained for x^^) at w = ± 1/4 the indicial exponents 
with a multiplicity of five (0^). While the diff-Pade analysis can only confirm this 
multiplicity of five as a lower bound it is the correct one at w = 1/4 as shown by our 
exact (conjectured) equation (51). 

From the linear combination of the long series for x''^'') X^^^ ^iiid x, we were able 
to make some serious progress on two important questions for a deeper physical and 
mathematical understanding of the full susceptibility x- 

Firstly, we have, finally, resolved the issue of the power/log behaviour of each 
X^"^ at the singular points versus the behaviour of the full x- The diff-Pade analysis 
performed for the full x shows the non-occurrence, at w = ±1/4, of the logarithmic 
singularities corresponding to the individual x^"^ , but shows with good accuracy the 
indicial exponents known to occur in x- We presented a model showing the mechanism 
of the resummation of the infinite number of logarithmic singularities that proliferate 
in the holonomic n-fold integrals x^"' building the known divergence of x o,t scaling. 

Secondly, we have proved by a Landau analysis that no non-Nickelian singularities 
of can lie on the |s| = 1 boundary of the principal disc. Since we do not 
understand the origin of the additional ODE singularities at w = 1/2 for x^^^ and 
vP' = 1/8 for x^^-* we cannot make the same analytic claim. Instead we have verified 
this numerically by fast-Fourier transform on appropriately modified series. The FFT 
results are surprisingly clear-cut, and display peaks at the precise points of Nickelian 
singularities on the unit circle. The possibility of cancellation becomes more and more 
remote, up to a point where we are able to confirm the existence of a natural boundary 
for the full susceptibility. 

In conclusion, we can say that, with this kind of experimental mathematics based 
on "extreme" massive computer calculations, we are exploring a new kind of "modular" 
lattice statistical mechanics, getting results that were unthinkable before. Indeed, we 
are now far along the road to a complete synthesis and understanding of the Ising 
model susceptibility. However, these experimental mathematical ideas for studying 
Ising model integrals are much more widely applicable, as they can be applied to any 
n-fold integral corresponding to a Feynman diagram. In such cases we are likely to be 
restricted only to the experimental aspects, leaving open the challenge for a subsequent 
theoretical understanding, such as that which we have been able to provide in the case 
of the Ising model. 



Extreme series 



36 



Acknowledgments 

We would like to thank B. McCoy for many illuminating comments. We would like 
to thank A. Bostan for helping to find the x^^^ ODE, modulo a prime, with his 
MAGMA program. We thank J. Dcthridge for some optimization of earlier C++ 
programs. The calculations presented in this paper would not have been possible 
without a generous grant of computer time on the server cluster of the Australian 
Partnership for Advanced Computing (APAC). Wc also gratefully acknowledge use of 
the computational resources of the Victorian Partnership for Advanced Computing 
(VPAC). One of us JMM would Uke to thank MASCOS, Melbourne for hospitality 
where part of this work was initiated and completed. This work is partially supported 
by a PICS/CRNS grant. This work has been performed without any ANR or ERC 
support. IJ and AJG gratefully acknowledge financial support from the Australian 
Research Council, and the hospitality of LPTMC, Universite de Paris 6 where much 
of this work was carried out. 

E-mail or WWW retrieval of series 

The series for the various generating functions studied in this paper can be obtained 

via e-mail by sending a request to I. .Jcnscn@ms.unimclb.cdu.au or via the world wide 
web on the URL http:/ /www.ms.unimelb.edu.au/^iwan/ by following the instructions. 

Appendix A. Singularity exponent sum-rules (Fuchs' relations) 

Suppose the order M and degree £) linear ODE, Lmd{S{x)) = 0, with Lmd given by 
(26), has singularities at x = Xi with multiplicities qi. Besides the true singularities 
of the generating function ^(a;), the x = Xi may be apparent singularities when 
M = Mq, the minimum order, or they may be M dependent spurious singularities 
when M > Mq. In any case the ODE is necessarily of the form 

Pd{x)-L,-{x--^)''-\S{x)) + --- = 0, (A.l) 

-^^a: = a; • ^ - So + ^{Qi ■ M - ^Qi {qi + 1) - Si)x/{x - a;,), 

i 

with Xi the roots of the degree D head polynomial Po (x) . One can explicitly verify 
that So is the coefficient of p^^^ in the indicial equation of degree M at a; = 0, 

- sop^"^ H + constant = 0, (A.2) 

and thus is the sum of M exponents at a; = 0. Similarly, Si is the coefficient of 
in the indicial equation of degree at x = Xi, 

- Si -p^'"^ + • • • + constant = 0, (A.3) 

and is the sum of the qi singularity exponents at x = Xi. 
If we set X = l/y, the transformed ODE (A.l) is 

PD{l/y) -Ly-iy. ±)^-\S{l/y)) + • • • = 0, (A.4) 

-^y = y • ^ + So - Xl*^* ■ M - ^Qi {Qi + 1) - Si)/{1 - yxi), 
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and we can identify — sq + {it' ^ — \ Qi {li + 1) — Sj) with s^, the coefBcient in 
of p^~^ in the indicial equation of degree M at y = 0, 

P*^ - Soo • P^"^ + • • • + constant = (A.5) 

and this coefficient is thus the sum of the M exponents at y = 0. This exponent 
equivalence is the sum-rule 

So + +Y^{si + ^qi{qi + l)) = M-^qi = M- D, (A.6) 

i, all i 

and is the starting point for our specialisations below. It is easy to show that (A.6) is 
equivalent to the usual Fuchs' relations [28, 29] which arc sum-rules on all exponents. 
The advantage of (A.6) is that it makes explicit the role of those singularity exponents 
that arise as solutions of the indicial equations of (typically) much smaller degree 
dictated by the multiplicity of the head polynomial zeros of the linear ODE. 

We now assume that the true singularities and their associated exponents have 
been determined. If M > Mq, then (A.6) is a constraint on the spurious singularity 
exponents and is of little interest. On the other hand if M = Mq, then there 
are no spurious singularities and the either true singularities or apparent 

singularities. Furthermore, in many cases of interest as we found in our analysis of x^^^ 
and x^^"* J the multiplicities of the apparent singularities at X — Xi circ all Qi = 1 and 
the exponents are all Sj = Mq. Then the contribution to ^^(sj -l-gj ((j'j-t-l)/2) coming 
from the apparent singularities is simply ^^^(si + l) = X^i(-^o + l) = Dapp- {Mo + 'i-) 
with Dapp = D —Dq and Dq the minimum possible degree of the ODE. The sum-rule 
(A.6) can now be rewritten as 

So + Soo + X) + ^ * • + = (^•^) 

i, true 

Mo-D - Dapp ■ (Mo + 1) = Mo • £>o - Dapp, 

which is an explicit formula for Dapp in terms of the true singularities only. To 
emphasize this point we note that for any physical problem the true singularities are 
determined and encoded in the generating function S{x). By deciding to represent 
this information as a linear ODE of minimum order of the form (26) we are forced to 
specify a total of (Mq + l)(£'o +Dapp + 1) coefficients with Dapp given by (A. 7). For 
the S = x^^^ example treated in the text, sq = 192 from Table 4 and Soo = 147, 
which is 56 from Tabic 4 plus the 14 term sum 0+1 + • • • +13 = 91 from the additional 
regular solution exponents inferred to be in the degree Mq = 33 indicial equation. 
The remaining true singularity exponent sum in (A. 7) is 653 from Table 4. On solving 
(A. 7) for Dapp we get the value Dapp = 1384 given in the text. The number of ODE 
coefficients is (Mo + l)(£'o + Dapp + 1) = (33 + 1)(72 + 1384 + 1) = 49538. 

There are a number of related observations that are significant. Firstly, the 
presence of apparent singularities in the minimum order ODE implies that there are 
constraint conditions [30], namely Mq — 1 conditions for each singularity in addition to 
the observed exponent value Si = Mq. Thus Dapp-Mo coefficients in total are fixed by 
the constraints and this means that in principle we need only {Mq + 1) (-Do + 1) + Dapp 
series terms to determine the minimum order ODE. In the x*-^^ example this is 
(33 + 1) (72 + 1) + 1384 = 3866 which is much smaller than the > 7000 terms 
for any x^^^ ODE given in Table 1. Unfortunately, we know of no practical way to 
implement the constraints as they are in general non-linear. Specifically, let us write 
the minimum order ODE operator as /^(a;- ^)^° with fm= for m > Mq 
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and the head polynomial /o factored as /o = Ptme " Papp^ thus clearly separating 
the true and apparent singularities. The apparent singularity constraint conditions 
are then the statement that each 

Pm ~ [fm ' {fm ^ ^ ' ) fm—1 ' (/m+1 ~H ' Pcipp^ 

m= 1, 2, • • • , Mo, (A.8) 

is a polynomial of degree 2-Dq +Dapp. It is the exact division and reduction in degree 
from 2 (Do + Dapp) that implies the existence of Dapp conditions in the numerator 
of the right hand side of (A.8) for each individual m. 

Secondly, although the number of coefficients needed in an ODE can be 
dramatically reduced by moving away from minimum orderjj, this is offset by a 
dramatic increase in the size of the integer coefficients specifying the ODE. This 
observation is based on our experience with a 6 x*-^' — X^^^ analysis. The minimum 
order ODE has Mo = 6 and degree = 40. With the normalization choice omo = 1 
in (26), the remaining coefficients arc integers and can be found by a mod prime and 
Chinese remainder theorem analysis using five primes. In contrast, with the same 
normalization, the "best" choice ODE from Table 2 with M = 10, D = VJ has as 
its remaining coefficients rational fractions with numerator and denominator integers 
each typically 160 digits in length. To find these requires a supplementary continued 
fraction calculation starting from Chinese remainder theorem residues of 320 digits and 
requires about 70 primes. Similar results hold for other non-minimum order ODE's 
for 6x('^) — x*^^) and we believe analogous disparities are likely in any 2x^^^ — X^^^ 
analysis. But we need to emphasize that the utility of the "best" mod prime ODE as 
a recursion relation device [24] remains. It enables us, for any given prime, to extend 
the shortest possible generating function series to the length necessary to find the 
minimum order ODE. 



Appendix B. On the Landau singularities 

The discussion of the Landau singularities of Ising like integrals in [14] is based 
on analogues of the (2n — 2)-dimensional integrals for x^"^ given by Wu et al. [2]. 
Here we derive and extend those results based on analogs of the (n — l)-dimensional 
integrals (5). The calculations are complementary, each having certain advantages 
and disadvantages. 

In the forthcoming technical discussion we will often use the following definitions 
supplementing those in (5-9). The square root factor appearing in a;, and we denote 
by: 

/i = ^ • V(l-2^«cos(0,))'-4«;2 = l/y,. (B.l) 

We also define 

cos(Ci) = l/{2w) - cos((/.i), sin(O) = i • /i (B.2) 

so that we can write: 

Xi = exp(iCi). (B.3) 

The formulae (B.1-B.3) are understood to apply for s or w; small and elsewhere 
by analytic continuation. Furthermore we will take it to be understood that by the 

(I In the x^^^ example we can move from the 49538 coefficients given above for M = Mq = 33 to 
(M + 1)(D + 1) = 7410 given in Table 1 for M = 56. 
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integral x^"^ we mean (5) and various analytic continuations of the series analytic at 
w — 0. This means wc do not distinguish between "integral" and that "particular 
solution" of the linear ODE that agrees with (5) for small w (or s) but whose domain 
is not restricted in any wayf. This is a great simplification for our discussion below 
since the value of x^"''(^c) reached by analytic continuation in general depends on 
the path chosen for w between and Wc and this can lead to very complicated 
topological considerations. We do not address any of that here but note that some of 
the differences in x^"^ on different branches might arise because analytic continuation 
requires the displacement of branch cuts which result in changes in the signs of fi and 
Q in (B.1-B.3). The results we describe below allow for all possible sign changes and 
thus the singularity list is the complete list covering all these local signs. 

This must be borne in mind when, for example, we say x^"^ is singular at certain 
points |s| < 1 since it is certainly the case that the radius of convergence of x^"^ is 
\s\ = 1. The latter observation follows trivially from the fact that if the phases in (5) 
are real then for |s| < 1 one finds that fi cannot vanish and also \xi\ < \s\. When 
we want to restrict the domain on which the integral is defined^ to the region |s| < 1 
(and on which for \s\ < 1 the series applies), wc will explicitly indicate this by making 
reference to the "principal disc" or the "integral on the principal disc" . 

The points s = ±1, ±i or equivalently = ±1/4, oo are fairly obvious singularities 
of X*'"'' and to simplify the analysis wc explicitly exclude these points. A consequence 
of this exclusion is that (B.2) implies that if sin(^j) = sin((/)j) = then w = ±1/4 or 
oo. Thus when we exclude these singular points we are preventing sin(Cj) and sm{(pi) 
from vanishing simultaneously, a result we will use on a number of occasions below. 

The discussion of the Landau singularities of an integral like (5) is somewhat 
simplified by the fact that the integration is over a unit cell of a periodic function 
and thus there arc no end-point singularities. In fact it is best to think of the 
integrations as closed contour integrations in Zi = exp(i0j) that can be arbitrarily 
deformed away from the unit circle provided no integrand singularities are crossed 
while the constraint H -^i — 1 maintained. These singularities in (5) are at /,; = 
for all i, XiXj = 1 for all i ^ j and H^^i = l- They are not all independent since, 
for example, fi = implies Ci = or tt and Xi = ±1. Thus the vanishing of fi 
can overlap with the Xi product singularities. However it is important to note that 
fi = 0, which is a condition on a single (pi, is very different from H^^i = l) which 
is a relationship between all ^j. If /, = cannot be avoided by the Zi contour for 
some particular w, then the integral x''"'' is singular at that w. On the other hand 
]^ = 1 can lead to a singularity of the integral only if Y[ Xi is also stationary with 
respect to variation in all cpi subject to the phase constraint in (9). Both = and 
]^ a;^ = 1 arc examples of pinch singularities, but to distinguish the simpler fi = 
case we will refer to it exclusively as a pinch singularity in the discussion below. The 
n = 1 case or any similar situation in which a non-trivial stationary condition must 
also be satisfied we will call a van Hove singularity in recognition of his analysis [31] 
that predates that of Landau [32] by several years. Another singular integral situation 
arises when, say, xiX2 = I and H = 1 are simultaneously satisfied but rather than 

t Thus completing the generaUzations already in print. In [10] "integral" was understood to mean 
a single-valued function in the cut w plane. The definition of "integral" was extended in [15] to be a 
single valued function in a cut s plane and hence double-valued in w. Finally, in a Landau singularity 
analysis as in [14] one no longer attempts to specify on which local Riemann sheet any particular 
singularity occurs. 
^ And is univalued! 
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X1X2 and Yl Xi being separately stationary the normals to these two hypersurfaces are 
parallel ([33], p. 48). In this case the integration variables are trapped between two 
distinct hypersurfaces that touch tangentially. 

Both to simplify the discussion and because this includes the most important 
situations, we start with the Landau problem of only ]^ a;, = 1 and fi = as 
singularities in the integrand of x^"^- Since Xi = exp(ii^i) this product constraint 
together with the phase constraint in (9) results in the symmetric pair 

^ = mod 27r, ^ = mod 27r. (B.4) 

The phase constraint we handle directly by taking 4>i,i = 1, • • • , n — 1, as independent 

so that Ci = C,i{4'i), « = 1, 1, and Cn = C,n{4>n) = Cni^nk - 4>i 0„-i). 

The requirement that be stationary with respect to phase variation is that the 
derivative combinations C,[ — ('^ vanish where i = 1, ■ ■ ■ , n — 1 and Ci = dQ/dcpi = 
— sin{(f)i)/ sin(Ci). The last equality follows from the definition (B.2). The derivative 
conditions can be rewritten as 

sin(Cj) • sin((/)j) = sin(Cj) • sin(4), i^j (B.5) 

and our derivation requires that (B.5) be subject to the restriction that no sin((^j) 
vanishes. However this restriction can be dropped because we are looking only for 
solutions for which sin(c6i) and sin(Ci) arc not simultaneously zero and the case that 
(B.5) yields all sin((^j) = is nothing but the pinch singularity condition we must 
investigate also. We can rewrite (B.5) for each i,j combination as the pair 

(cos((^i) — cos(^j)) • {4w — cos{^i) — cos{^j) + 4:W ■ cos(0i) cos(0j)) = 0, 

(cos(Ci) - cos(Cj)) • (4w - cos(Ci) - cos(Cj) + 4w • cos(Ci) cos(Cj)) = (B.6) 

by squaring, rewriting the sine functions in terms of cosines, and utilizing the definition 
(B.2) which we reproduce here in symmetric form 

cos((/>i) +cos(Ci) = l/(2'u;), i = 1 . . . n. (B.7) 

Note that (B.6) and (B.7) allow all possible sign combinations ±(/)j, ±(j for a given 
(pi, Q pair and one must in all cases check that (B.5) is also satisfied. That there is a 
remaining sign degeneracy allowed by (B.5) is a consequence of our decision to define 
our integral to include all possible sign combinations in the local Riemann sheet. If 
the domain of w were to be restricted to a particular Riemann sheet then additional 
analysis would be required to determine the uniquely signed solution. 

The equations (B.4-B.7) are the Landau conditions for our reduced problem of 
Ha^i = 1 and fi = as the only singularities. We now consider specific situations, 
Case 1 to Case 5. 

Appendix B.l. Case 1: the phases 4>i and Q equal or w 

The simplest case is that of phases (pi and Q equal or tt. The stationary constraint 

(B.5) is trivially satisfied and (B.7) implies singularities at w = ±1/4 or 00. These are 
the points we excluded from the analysis and we again note that sin((/)j) and sin(^i) 
cannot simultaneously vanish. 

Appendix B.2. Case 2: all (pi equal and all Q equal 

Equally simple is the case of all (pi equal and all d equal. Again the stationary 
constraint (B.5) or (B.6) is satisfied automatically and the phase constraints (B.4) are 
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satisfied with 

= 2'Kk/n, (i = 27rm/n, < k, m < n (B-8) 

and with and Q not both or tt so as to exclude Case 1. The associated singularities 
are at 

l/{2wk,n) = Skm +l/sfem = cos(27rfc/n) + cos(27rm/n). (B.9) 

Because the stationary constraint for these singularities is automatic the full Landau 
formalism is not necessary and (B.8), (B.9) could have been guessed just as they were 
in [3]. Note that |sfe„i| = 1 and the calculation in [3] shows them to be principal 
disc singularities. In the following we will designate them as Nickelian singularities 
to distinguish them from other van Hove singularities. The singularities Wkm can be 
given as the roots of polynomials in w with integer coefficients. 

Remark 1 A new situation not considered in [3] arises out of the Landau 
formalism. For a given <;ij, Q combination, the constraints (B.5-B.7) allow 4)j, Q 
to be — Ci ill addition to the (pi, Q we have considered. For every such sign 

reversal there is one pairwise cancellation in the constraint sums (B.4) so that we 
should add to (B.8), (B.9) new singularity conditions obtained by the replacements 
n^n — 2, n^n — A, etc. We will however leave (B.8), (B.9) unchanged and 
if there is a possibility of confusion, refer explicitly to (B.8), (B.9), for which there 
has been no pairwise cancellation, as the "irreducible" conditions. The n — > n — 2m 
replacement singularities will be treated separately as Case 5 below as this situation 
arises numerous times. 

Remark 2 A very important observation concerning the n — > n—2m replacement 
singularities is that none of them are principal disc singularities. To sec this note that 
(B.l), (B.2) defines Q as an even function of (pj for \s\ < 1 on the principal disc and 
by continuity to the limiting case |s| = 1 as well. On the other hand, for fixed ^i, Q, 
(B.5) requires C,j to be an odd function of (pj. Thus while (B.5) allows 4>j, Q to be —(pi, 
~C,i in addition to the (pi, this singularity cannot be on the principal disc but rather 
must be on those other "Riemann sheets" on which the fj square root function has the 
opposite sign. The numerical evidence is consistent with this result. It was already 
shown in [3] that x^^^ principal disc singularities were not present on the principal 
disc of x'^' . 

Appendix B.3. Case 3: all sin(Ci) = 0, i.e. all fi = 0, so that Ci = or it 

At the next level of complexity consider the possibility that all sin(^i) = 0, i.e. all 
fi = 0, so that Ci = or tt. The constraint (B.5), which in this case specifies a pinch 
singularity, is trivially satisfied. We find from (B.7) that there arc only two possible 
values for cos((p i), namely cos(0(+)) = l/(2«;) + 1 and cos(0(-)) = l/{2w) - 1. The 
former is associated with d = n and if we demand that the singularity condition 
in (B.4) also be satisfied then we must have an even number of these terms. We 
implement the phase constraint (B.4) on (pi as follows. Since there are only two 
possible valuesi we can write k(p^~^^ + (n — k)(p^~^ = mod 27r and from the 
preceding remarks, k is even. Equivalcntly, exp(ifc(/)'^"^^ + i(n — k);^*^^^) = 1 or 
exp(iA:0*^+)) = exp(— i(n — k)(p^~'>). Now add to this last equation the reciprocal and 
obtain cos(fc(/)("'")) = cos((n — k)(p^~^). This form is convenient because each cos(m0) 

X Exactly as in Case 2 we consider only the "irreducible" case in which a given value of cos{(j>i) 
defines a uniquely signed 4>i and hence there are no pairwise cancellations in the constraint ^ 4>i. 
We will remark further on the general situation, considered as Case 5, later. 
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can be expressed simply in terms of cos{<p). On using the definition of Chebyshev 
polynomials we get 



as the defining equation(s) for the singularities. Note that we have excluded fc = or 
n since all (pi are equal and (B.IO) becomes a special case of (B.8), (B.9). 

Remark 3 The singularity conditions (B.IO) do not allow for singularities on 
\s\ = 1, as we now show. The condition \s\ = 1 is the condition w real and 
—2 < l/(2w;) < 2. Consider first the interval < l/{2w) < 2. Then the argument of 
Tk in (B.IO) lies between 1 and 3 and \Tk\ > 1, while the argument of T„_fe lies between 
— 1 and 1 and jl^n.^l < f . Thus (B.IO) cannot be satisfied. A similar argument applies 
on the interval —2 < l/{2w) < 0. Since we excluded w = ±1/4 and oo at the outset, 
the proof is complete. 

Remark 4 The same integral singularity conditions (B.IO), but with odd index k 
allowed as well, were given in [16] based on the integrand singularity condition fi = 
for all i irrespective of Jl ■ We see here the partial overlap between the all ft = 
and n = 1 conditions, but because fi= is the more general condition it would 
appear that odd k should be included in (B.IO). In fact (B.IO) is correct because of 
the very special structure of the Ising integrals (5). This will be proved below. But 
first note that the integrals evaluated in [13], where the integrand in (5) was replaced 
by the simpler pole product (Hyi) j 9^11 show that both odd and even k are required 
in that case. That is to say, for these integrals the pinch singularity conditions are 
both necessary and sufficient and the same should apply also to the Ising x^"^- The 
only way contours that are trapped in the process of analytic continuation in the toy 
integrals could escape being trapped in x^"^ is for at least one fi factor to become non- 
singular. The mechanism for this to happen is that in the vicinity of a potential pinch 
an fi appears not in a Laurent series but rather in a Taylor series in powers of (fi)"^. 
Given the complexity of the Ising x^"-* integrand this would seem rather miraculous but 
in fact happens as we now show. Not too surprisingly, this point was entirely missed 
in [16]. 

The argument below is for a pinch singularity defined by /j = for all i. We 
begin by noting that at a point where fi is small, Q is either small or near tt. For 
purposes of the present argument only, we replace those near n hy d — w and 
incorporate this change consistently by the simultaneous replacement Xi — ^ —Xi. The 
factor (1 + n / (1 " n ^i) '^^ X^"^ integrand can then be written as cot (^ Ci/2) 
if an even number of such replacements were made or tan Ci/2) if an odd number 
of such replacements were made. In the first case, only a Laurent expansion is possible 
and the pinch singularity is qualitatively like that in the toy integral with integrand 
(n^/i)^- Since that integral was observed to have the singularities (B.IO) the same 
must be true here and we conclude that the necessary conditions (B.IO) are also 
sufficient for an integral singularity. 

However, in the odd replacement case a Taylor expansion in odd powers of ^ Q 
exists. This will yield a sum of terms of the form Y[ (C:)''' with with '^Pi odd. Since 
fi = i - sin(^j) implies Ci is an odd function of fi the expansion could equally well be 
written as a sum of {fiY' with again ^Pi odd. Because J^Pi odd there is at 
least one pi which is odd and of course this pi > 0. Now recall that the x*^"-* integrand 
also contains the factor J^j/j = n(l//i) ^^'^ t^is shifts every pi down by one. In 
particular the original odd pi has now become even and is still non-negative so that 



Tk{l/{2Wk) + 1) = Tn-k 




(B.IO) 



< fc < n. 
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the integration over the analytically continued <pi encounters no singularity, i.e. no 
pinch. This completes the argument for those proxy Ising integrals that do not contain 

the (G("))^ factor. 

To show that the presence of (G^"^)^ in the integrand does not change the 
argument for an absence of singularities we first note that G^"^ is a product of terms 
each of which is of the form 

2 sin((.^, - cPj)/2) ■ ^ff^ (B.ll) 

The XiXj dependent factor is proportional to either sin((Ci+Cj)/2) or 1/ s\n{{C,i+Qj) /2), 
after making the replacements described in the penultimate paragraph, so that now C,i 
and C,j are understood to always be near zero. The sine function in the numerator of 
(B.ll) is crucial in canceling a potential zero in the denominator, but has a complicated 
representation in terms of Q and Cj. For purposes of a singularity existence analysis it 
is legitimate to replace it by any more convenient analytic and asymptotically linear§ 
function in the neighbourhood of Q = Q =Q and we choose instead of sin((</)i — 4'j)/2) 
the difference (fj)^ — {fif = sin^(Ci) — sin^(Cj). Then (B.ll) is proportional to one of 

(sin2Ci-sin2g)sin((Ci + C,)/2), 

(sin^ Ci - sin^ 0)/ sin((C, + 0)/2) = (B.12) 
2 {sinQ - sinQ) ■ cos((C< - 

Both of the expressions in (B.12) have a Taylor expansion with terms of the form 
{CiT* [QjY' with Pi +pj odd. An equivalent expansion is in terms (/i)^* {fjY^ with 
Pi + Pj again odd. The square of this expansion will be similar but with pi + pj now 

even. Finally the expansion of (G*-"-*) must be a sum of terms of the form Y\ [fif^ 

with YliPi 6ven and multiplying the expansion obtained in the absence (G^"')^ by 
terms Y\ {fiY* for which J^Pi even cannot change the argument for the absence of 
a pinch. 

To summarize this case, the pinch conditions give (B.IO) as necessary. The 
evidence from the toy model study [13] is almost certainly a proof that (B.IO) is 
also sufficient for these singularities of x*^"^ . 

Remark 5 A situation closely related to sin(^j) = 0, for all i, as treated above 
is sm{(l)i) = 0, for all i. This leads to van Hove singularities in which the stationary 
conditions (B.5) and both conditions (B.4) must be satisfied. Each (pi is either or tt. 
The phase constraint then requires an even number of 0i = tt values. On imposing 
the singularity constraint in (B.4) the same solution (B.IO) is obtained. 

Remark 6 This equality of solutions for all sin(Ci) = and all s'm{(pi) = 
is expected and illustrates a very important aspect of the Landau analysis of the 
integral representation (5) for x*'"^. The original representation for x^"^ given in [2] 
is an integral over two sets of r?, — 1 independent phases with an integrand that is 
symmetric under the interchange of these sets||. Integrating out one set to arrive 
at (5) has obviously broken this symmetry, but a vestige of it remains in the two 
conditions (B.4). The 4>i are the remaining dummy phase variables of integration 
while the d = Ci{4'i) ^ii'c the specific values that result from evaluating the other set 

§ Linear in the deviation of (pi from the singular point. This means a Hnear function in tiie quadratic 

II The invariance under phase variable interchange is not always obvious in a particular formula. See 
the remarks in Wu et al. [2] following their equation (4.87) on this point. 
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of phases at the Wu et al. [2] integral singular points. Finding singularities of the 

integral such as (B.IO) amounts to completing the integration process and while the 
steps associated with dealing with the "integrated" Q versus the "unintegrated" 
are obviously different, once the (pi values have been fixed by the Landau conditions 
the symmetry is restored and the final answers one obtains must be the same. That is, 
for every combination of the 4>j leading to an integral singularity there is another 
set (f)i, Q obtained by <-> interchange that leads to the same singularity and we 
can always choose the combination that requires the least computational effort. 

Remark 7 The above comparison of alternatives also illustrates that there is no 
fundamental distinction between the pinch and van Hove singularities. Rather it is just 
a distinction in the ordering of the steps of the calculation. This is particularly relevant 
when comparing the calculations here to those in [13]. There the Landau analysis is 
done at the (2n — 2) integral stage and there is no analogue of the pinch singularity in 
the restricted sense we have defined here. Instead, every singularity evaluation requires 
a supplementary stationary condition evaluation and thus all singularities are of van 
Hove type. 

This concludes our treatment of the singularities associated with the pinch 
condition sin(^j) = and the corresponding special van Hove case sin(0j) = 0. 
This list is exhaustive as there can be no integral singularities for the mixed case that 
sin(Cfc) = 0, fc = 1, • • • , m, and sin(Cj) ^ 0, j = m + 1, • • • , n, with I < m < n. 
To see this note that while the (k is a singular, i.e. square root, function of (pk in the 
neighbourhood of sin((^fe) = 0, the (pk is an analytic function of Cfe- Specifically, (pk is 
asymptotically a constant with an added quadratic dependence on the deviation 5(^k 
of (^k from or tt. If, for each fc, 1 < A; < m, we replace / d(pk/ sin(^/5) ~ / d<pk/fk by 
— / d(k/ sm{(pk) then the only integrand singularity condition remaining is ^ (^j = 
mod 27r. Furthermore, because we are explicitly demanding sin(Cn) 7^ 0, Cn{<Pn) is 
an analytic function in the neighbourhood of the integrand singularity. As before 
we eliminate (pn using the constraint = mod 27r except that here we also 

express (p^ as (pkiCk) for 1 < k < m. The leading dependence of Cn on (k is the 
quadratic dependence 6C,^ and this implies that no variation such as d^Q/d(k = 
dQk/dQk — d(n/dCk = 1 — constant (5CfeC can ever vanishf. That is, no integral 
singularity is possible in the mixed case. 

Appendix B. 4. Case 4: neither sm{(j)i) nor sm{Q) vanish 

In all of the above cases the Landau stationary condition was either automatic because 
of symmetry or did not need to be invoked because the same singularities arose from 
the fi=0 for all i condition. References [13, 14] showed how other singularities can 
arise out of the Landau rules for some toy integrals and we give in this section the first 
and simplest example of a case where the Landau stationary conditions are non-trivial 
for the Ising x''"^- Consider the situationi that in (B.4-B.7) k of the cpi, Q take on 
the values (j)a, Ca and the remaining {n — k) are (pb, (1,. The phase constraints can be 
written as polynomial relations in cosines as in (B.IO) but with the difference that the 
individual terms like cos{(j)a) are not explicitly given as functions of w but are to be 
determined self-consistently together with w. In summary we must solve the analogue 

t That Cn is finite is crucial. This fails in the previous pinch situation where all sin(^i) = 0. 

X As in Cases 2 and 3 we restrict our attention to the "irreducible" situation in which there is no 

cancellation between pairs of phases differing only in sign. The general situation is treated as Case 

5. 
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of (B.IO) which is the pair 

Tk{cOS(l)a) = Tn-k{cOS(j)b), 

TkicosCa) = T„_fe(cosC6), < fc < n (B.13) 

together with the stationary constraints (B.6) 

4w — cos (l)a — cos 06 + ■ COS (/Jjj COS (ia = 0, 

Aw — COS (a — cos (b + Aw ■ COS Cb COS (a = 0, (^-14) 

and the definitions (B.7) 

cos (j)a +C0sCa = l/(2w), COS (pb + COS = l/(2w). (B.15) 

We can reduce these equations to somewhat simpler form by using (B.14) to 
ehminate <pb and ^5- The result is a triplet of equations 

Tk{cOS(()a) = Tn-k{{Aw - COS (t)a)/ {I- Aw COS (pa)), 

7fe(cosCa) = Tn-k{i4:W - cosCo)/(l - 4wcosCo)), < k < n, 

cos (t>a+ COS (a = 1 / {2w) (B.16) 

defining the van Hove singularities when neither sin(0j) nor sin(Cj) vanish. 

To produce the singularity polynomials, we use the last condition in (B.16) to 
eliminate cos(^a). Then, the elimination of cos{(f>a) from the first two polynomial 
conditions in (B.16) gives the roots that are inserted into (B.16) to determine which 
are actually solutions. We must also go back and verify that (B.5) is satisfied as our 
candidate polynomials are all based on the squared form of (B.5) as expressed in (B.6), 
(B.7). 

Our explicit procedure for eliminating the phase variables in (B.16) is as follows. 
The first condition in (B.16) can be expressed as a polynomial equation pi(cos(0o)) = 
0. The second condition can be reduced to another P2{cos{4>a)) = if we use 
the last condition to eliminate cos(Ca)- We are not interested in the solutions 
cos{(f>a) = ±1 as these are a subset of those already found. We therefore take it 
as given that pi and P2 are reduced polynomials which do not contain the factors 
cos{4>a) = ±1. They are both of degree m = n — 2, n — 1, n, depending on 
whether k and n — k are even-even, odd-even or odd-odd combinations. We now 
eliminate cos{(j)a) by the following iterative process. We first reduce the algebraic 
complexity of the polynomials pi and p2 by dividing them at each stage by their 
greatest common divisor. Assume this has been done and they are now written as 
Pi = ai + ■ ■ ■ + bi cos™ (0a), P2 = a2 + • • • + 62 cos"^{(j)a)- We then generate two new 
polynomials qi = {A2P1 — A1P2)/ cos(0a) and 52 = B2P1 — BxP2 where Ai and A2 
arc just fli and 02 divided by their greatest common divisor and similarly for Bi and 
i?2- By construction gi and q2 are of degree one less than px and p2 and after division 
by their greatest common divisor are the p\ and P2 of the next iteration. The process 
stops when we reach degree 0. The greatest common divisors are polynomials in w 
whose zeros are candidates for solutions of the original equations. 

Remark 8 We find empirically that for given n and A; in (B.16) the degree of the 
w polynomials giving the Landau singularities is bounded by n{n — 2)/A — m.{m,— l)/2, 
m = min(A:, n—k). There is evidence from [14] based on toy integrals that the predicted 
singularities are seen in the ODE so that just as in Case 3 we conjecture that our 
polynomial solutions for Case 4 are both necessary and sufficient. For n — 2N + 1, 
the number of s plane singularities is N{N + 3), N{2,N + 1) and bN{N'^ - l)/3 for 
Cases 2-4 respectively. The corresponding number of s plane singularities for n = 2N 
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is 2[{N - 1){N + 3)/2], 6[{N - 1)2/2] and 2[{N - 2){2N - 3)(5A^ + 4)/12] barring 
accidental degeneracies§. 

Remark 9 As continuation of our empirical observations we note that of the 
Case 4 singularities, N{N - 1) are on \s\ = 1 iov n = 2N + 1 and 2\{N - 2f /2] 
for n = 2N . However, in none of these cases is the singularity on the principal 
disc of x^"' and the general proof that no Case 4, \s\ = 1, singularity lies on the 
principal disc is as follows. Observe first that a singularity with \s\ = 1 on the 
principal disc requires that the integrand of x'^"'' be singular for real phases 0^. 
That is, both (pa and 0(, must be real. Equivalently cos(0o) and cos((/>(,) are real 
and of magnitude less than one. Furthermore, the cosines are related by (B.14) or 
cos((/)b) = (4u> — cos((/)a))/(l — 4wcos((^a))- If 4w is real and \4w\ > 1, which is the 
condition for \s\ = 1, then |cos(0o)| < 1 yields |cos((/)(,)| > 1 and vice versa. An 
important consequence of this proof is that there are no Case 4 singularities to cancel 
any of the Case 2 Nickelian singularities. 

Although we started the discussion of Case 4 with the remark that our calculation 
was just the first and simplest example of a non-trivial van Hove singularity, one 
observes that the singularity conditions (B.6), for example, require cos((/)j) — cos(c/)i) 
or cos((/>j) = (4w — cos(0i))/(l — 4wcos((/>j)) and the latter is an example of a one to 
one Mobius' mapping. Thus it is not possible, for any given w, to have more than two 
distinct values in the set cos{(f>i), i = 1, • • • , n. This in turn implies that as long as we 
consider f| a;j = 1 as the only singularity we have exhausted all possibilities for the 
values of cos{^i) in Cases 2-4 above. Of course this applies only to the "irreducible" 
situations and we turn now to the implications of the and ±(i sign degeneracy. 

Appendix B.5. Case 5: ±0i and zLQ sign degeneracy 

Once the "irreducible" solutions Cases 2-4 axe known, all other solutions can be 
obtained by recursion which we now describe. Suppose the combinations 5^"^ = 
|(0|"\Ci"^)) i = li ■ ■ ■ ) '^1 are the stationary singular points of the integrand of 
-^(")_ Then we can construct a stationary singular point of the x'""^^'' integrand 
as S'("+2) = |s'("),(</-j"\cj"'),(-0j"\-Cj"^)}, where j is any one of the n values 

in S^^\ Clearly by adding a pair with opposite signs wc have guaranteed that the 
phase constraints (B.4) are satisfied. Furthermore, all constraint equations involving 
only cosines remain unchanged and thus satisfied. Finally, the stationary condition 
sin(^i) sin(0j) = sin(0i) sin((^j) is potentially satisfied under the simultaneous sign 
change of and Q . Whether this occurs will depend on which Riemann sheet one is 
on. As an important example, and previously discussed in Case 2, on the principal disc 
for \s\ < 1 the definitions (B.l, B.2) show that sin(Cj) is an even function of 4>j and the 
new stationary condition generated by recursion cannot be satisfied. Specifically, the 
recursion mechanism explicitly excludes the possibility of addition to, or cancellation 
of, Case 2 "irreducible" singularities on the principal disc and leaves the argument 
that \s\ = 1 is a natural boundary of x secure. 

In conclusion, the Landau conditions for a singularity of x'"'' are also the 
conditionsjl for that singularity in ^^("+2^)^ fc > 1. Thus for a complete picture of the 

§ Af = 6 affords examples. Here the [fc,m] combinations [0,3] and [2,2] in cos{kn/N) + cos(m7r/Af) 
yield the same l/(2w) = 1. Other degenerate combinations are [0,4] and [2,3] and those obtained 
by the replacements fc^6 — fc, m— >6 — m. 

II As always, these conditions axe necessary but not sufficient. The numerical evidence from the 
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singularities of x^"^ it is sufficient to list the polynomials defining the "irreducible" 
singularities for the Cases 2-4. 

Appendix C. Singularities 

The following is a listing of the polynomials defining the Landau singularities of x^"^ 
for small n for the three Cases 2-4 discussed in Appendix B. Our notation for the 
polynomials from x^"^ is P{'-'°-'^'^ n) or P{^"'^'^ nn-k,k) with k defined in eqns. (B.IO) 
and (B.16). We do not include the physical singularities w = ±1/4 or the simplest 
unphysical singularities at w = or oo. We also do not give the (reducible) Case 5, 
which would just be the list of m > 0, polynomials for any given x^"''- An 

example of the notation we use in the text in this situation is P(^7/^3i^2) to indicate 
a Case 3, n — 2m = 3 contribution to x^^^ • 
Case 2. Circle singularities P(^n„,o) = -P(^n). 

P(2 3) = (l + 2w)(l-w), 
P(2 4) = (1-4^2), 

P(2 5) = (1 + w){l + 2w; - 4«;2)(1 - 2,w + w^), 
P(2 6) = (1 - Aw^){l - W)(l - w^), 

P(2 7) = (H- 2w; - 8u;^ - 8w^)(l + 2w-w^ - w^){l -bw + - w^), 

P(2 8) = (1 - 8u;^)(l - 2w^){l - 4«;2)(l - l^vi^ + 4«;^), 

P(2 9) = (1 + 2w)(l ~w){l- Uvj^ + 8w;3)(l - 6w; -|- 9^^ - w^) 

x(l - 2,w^ - w'^){l + Zw- V?), 
PC 10) = (1 - w'^){l - 5w^){l - 12w^ + 16w^) 

X (1 - 7w'^ + w'^){l - I5w^ + 25w'^). 
Case 3. fi = or cos<^+ = l/2w + 1, cos(/>_ = l/2w — 1. 
P(^3i,2) = (l+3w; + 4«;2), 
^(^42,2) = 1, 

P(3 53,2) = (1 - 7w + 5w;2 - 4w;3), 

P(3 5i,4) = (1 + 8w + 20w^ + 15^3 + 4w^), 

PC 64,2) - (1 - 25w2)(i _ + 16w/), 

PC 75,2) = (1 - Ww + '35w^ - 51w'' + 21^4 - 4w'^), 

P(3 73 4) = (1 + 7y, + 26^2 -I- 7w^ + Aw^), 

P(3 7^ g) = (1-1- 12u; -I- 54u;2 -|- 112u;3 + i05w^ -|- 35w^ + 4w^), 

P(3 86,2) = (1 - 16w;^ + 100u;^)(l - 20^^ -|- 16w^ - mw^), 

P(3 84,4) = (l + 2?«2). 

Case 4. /j ^ 0, cos (j)a and cos^^ distinct. 

P(4 32,l) = l, P(4 43,l) = l, P(4 42,2) = l, 

P(4 54,1) = (1 - w - 3«;2 4w^), 

singularities of the associated ODE is that they are also sufficient and each singularity occurs on 
some Riemann sheet of the integral. 



Extreme series 



48 



^(^53,2) = (l+4u; + 8w2), 
P(4 65^1) = (1 - lOw^ + 29w^), 

P('64,2) = l, 
P('63,3) = l, 

P(4 76,i) = (1 - 3w7 - 10?«2 + 35w^ + - 62w^ + 17w^ + 32w^ - 16m;®), 

P(4 75 2) = (1 + 8w + 15u;2 - 21^3 - eOw"^ + 16w^ + 96w^ + 64m;^), 

P(4 74 3) = (1 _ 4u, _ 16^2 - 48u;^ + 32u;* - 128«;^), 

P(4 87^1) = (1 - 26m;^ + 242u;^ - 960«;^ + 1685u;^ - 1138m;^°), 

P(4 86,2) = (l-10w;2 + 32«;4), 

85.3) = (1 - 30u;2 + 56u;^ - 1312w<^), 
PC 84,4) = I. 

Appendix D. The one term Fermionic toy model 

Here we consider the question of which integral singularities might arise from 
combinations of XiXj = 1, = and Yl^'i = 1- 

The analysis of x*-"'' with the full Fermionic factor (G^"))^ is complicated in large 
part because of the very many different Xj Xj = 1 combinations to consider. Thus we 
are motivated to look at simpler integrals intermediate between and x''"'\ but 

which lead to the same singularities seen in the ODE for x'-^'^ and x'-^K The simplest 
situation is where the full Fermionic factor is replaced, for n even, by (/ii2)^, and for 
n odd by (/i2)^, where [3]: 

/12 = ^(sin(/>i -sin(ji2)(l +a;ia;2)/(l -XiX2)- (D.l) 

The integrands in these cases are simple enough that long series can be derived with 
only slightly more effort than for ^^^^ ■ We find, for example, by a mod prime analysis 
of the ODE for this "one term" reduction of x^^^ , that the head polynomial has the 
(1 — 2w) factor that is in the x'""^' ODE and in addition has the factor (1 — 8w^). 
The latter also shows up in the "one-term" reduction of x^^\ where it might have 
been expected since that it is what wc found for the x^^^ ODE by diff-Padc analysis. 
Explicit solution of the ODE shows the associated singular functions have, as leading 
terms, (1 — 2iv)^^'^, ln(l — \/8w) and ln(l + \/8w) for the "one term" x^^^ and 
(1 — 8t(;^)ln(l — 8w^) for the "one term" x'-^K The observed singularities of the 
ODE make the simple "one term" reductions of x*'"' ideal integrals for a Landau 
singularity analysis and the details of this analysis is the content of the rest of this 
appendix. As in Appendix B, we allow for various^ analytic continuations of the 
integrands. Yet in spite of this we find that the Landau conditions are never satisfied 
at either 1 — 2w = or 1 — 8w'^ = 0. Since the Landau conditions arc necessary for 
the integral to be singular we must conclude that the linear ODE generated from the 
series have additional singularities not possessed by the integrals. 

The addition of the new singularity X1X2 = 1 requires an analysis that can be 
broken into two parts. The simplest is the determination of the integral singularities 

^ "Various" includes all possible sign combinations of the square root factors in the integrand and 
that means we are considering every possible local environment in the Landau analysis. The same 
local conditions can differ in the global behaviour of the contour distortions used in the integrals but 
this does not aifect the singularity conditions. In this sense our singularity search is exhaustive. 
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arising out of just X1X2 = 1 and /j = 0. The more involved investigation is for the 

combination of fi = with xi X2 = 1 and Y[ = 1 simultaneously satisfied and 
with the normals to the latter two hypersurfaces constrained to be parallelf. 

For the simple case of X1X2 = 1 and /« = for all i we first observe that the 
fi ~ Q pinch conditions lead to the Case 3 singularities (B.IO) and the only new 
investigation to be done is to determine whether the presence of different integrand 
factors changes the even k condition in (B.IO). The previous proof that the presence 
of (G^"') did not change the even k condition in (B.IO) relied on the fact that 
G*^"^ is a product of hij factors and thus was in essence also a proof for any single 
(/ii2)^. Essentially the same proof will apply to (/i2)^ provided we take for fij the 
original fij = 1/2 (sin(^,) — sin((?!)j)) (1 + XiXj)/{l — XiXj) (see (4) in [3]) because 
this form can be reduced to a Taylor expansion with terms (CiYHCj)^^ with pi + pj 
odd exactly as in the argument for hij following (B.ll). The cumulant reduced form 
fij = {sin{(j)i) — sm{(l)j))xiXj/{l — XiXj), which is more convenient for calculating 
the full ;^(2n+i)^ when expanded has terms of both odd and even parity and thus will 
also generate odd k singularities (B.IO). This has been observed and is an additional 
confirmation of the argument that it is the very special nature of the integrand of the 
Ising x^") that is responsible for the even k condition in (B.IO). It is to be understood 
in the following that wc will be using the original fij and thus that (B.IO) will not be 
supplemented with odd k terms. 

We must also deal with the van Hove singularities arising from X1X2 = 1 
with X1X2 stationary. The product X1X2 — 1 requires Ci = ~C2 and hence 
cos(Ci) = cos(C2)- This in turn, because of (B.7), gives cos(^i) = cos(^2) and 
(pi = ±(/)2. The stationary condition is sin((/)i) = sin((/)2) = 0, which combined 
with the preceding = ±02 allows as possible 0i, (j)2 combinations 0,0 and 7r,7r. 
The corresponding cos(Ci) and cos(C2) are both either l/{2w) — 1 or l/{2w) + 1 so 
that the equality cos(Ci) ~ 008(^2) required for X1X2 = 1 is automatic and gives no 
constraint on w. The absence of a first order stationary constraint requires that we 
go to second order so that in addition to ([ = — sin((/)i)/ sin((i) = we demand 
C" = - (cos((;Ai)sin2(Ci) + cos(Ci) sin2(0i)) / sin3(Ci) = and similarly for (2- The 
expression in braces factorizes into (cos((/ii) + cos(Ci)) (1 — cos(0i) cos(Ci)), which, 
given (B.7) and the possible values cos(0i) = ±1, reduces to l/'u;(l ± l/(4'u;)). Thus 
from = we get as the only possible Landau singularities w = ± 1/4 or 00, the 
Case 1 singularities we are not considering. 

We begin the analysis where X1X2 = 1 and H^^i = 1 are simultaneously satisfied 
and these hypersurfaces touch tangentially, by deducing the necessary constraint 
conditions. As before we treat (pi, i = 1, • • • , n — 1 as independent giving Q = Ci{(j)i), 
i = 1, • • • , n — 1, and C„ = Cn{(t>n) = Cn (27rfc — (pi — ■ ■■ <^„_i). Parallel normals 
requires ad {Ci + ■ ■ ■ + Cn) /d^i = /39(Ci + C,2)/d(pi, i = 1, • • • , n — 1, with both a and 
(3 non-zero. Explicitly, 

«-(ci -C) = /3-ci, «-(C2 -c;) = /3-C2, (D-2) 

a • (Ci - C) = 0, 3 < i < n (D.3) 

where again C,' = dC,/d(p = — sin((^)/sin(^). Conditions (D.2) imply either = 

t We need not consider the third possibiUty, X1X2 = 1 and f^i^i = 1, satisfied with 2:1x2 and 
Y\ Xi simultaneously stationary. The constraints when xiX2 and Y\ are treated as independent are 
either compatible with each other or not. If the constraints are incompatible solutions will have to 
be dropped, otherwise they can be kept. In either case there is no possibility of new solutions being 
generated. 
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or = (2 while conditions (D.3) reduce to the previous stationary conditions 
(B.5) except that i and j are restricted by i, j > 2. The new structure of the 
stationary conditions has imphcations for the pinch singularities also. Whereas 
in the discussion at the end of Case 3 we noted that pinch and van Hove type 
singularity conditions could not mix, here we have the possibility that with m = 2 in 
that discussion ad{Ci + --- + Cn) /dC,i = /35(Ci + C2)/3Ci and (Ci + • • • + Cn) /5C2 
= + C,2)/dC,2 can both be satisfied with a = (i. Thus it is possible to have 

sin(Ci) = sin(C2) = in conjunction with the conditions (D.3) that arise out of 
the remaining derivatives with respect to (/>j, i > 2. The reverse situation in which 
sin((^i) 7^ 0, sin((^2) 7^ and sin(^j) = 0, i > 2, is not possible because there is no 
analogue of the a = (] solution in this case. In summary, all possibilities covered by 
(D.2, D.3) plus the allowed pinch situations are given by 

sin((/)„) • (sin(Ci)sin((/)2) - sin(0i)sin(C2)) = 0, (D.4) 
sin(Ci) • sin(0j) = sm{(j)i) ■ sm{Q), (D.5) 
i^j>2, (sin(Ci)^0,i>2) 

Details for the various alternatives are given below with the major categories a) 
through d) being the distinctions allowed by (D.4). 

a) Pinch case: sin(^i) = sin(^2) = 0. 
Allowed combinations of Cii C2 satisfying sin((^i) = sin((^2) = and X1X2 = 1 
are 0,0 and 7r,7r. The corresponding cos((/)i) and cos((/)2) are either both cos((/)'^^)) = 
l/(2'u;) — 1 or both cos((/>^"'"^) = l/(2'u;) + l. = —(^2 then conditions (D.5) together 

with (B.4) are just the conditions for the singularities of x^"~^^ the absence of the 
Fermionic factor (G^"^) and have already been discussed. If (pi = <j)2 the sum on 
phases 0,, i > 2, is no longer 0, but ±2^i mod 2-i: with 0i either or 0^+). For 
the new singularities that arise we introduce further subdivisions depending on the 
form taken by (D.5). 

For i > 2, all equal and all (pi equal as the analogue of Case 2. The 
first new situation arises for x^^\ Here we set Q = 2tt/3 and solve T2{(p^'^^) = 
T3(l/(2w) + 1/2). For x^*^^ set Q = 7r/2 and solve T2((/)(±)) = r4(l/(2w)). These are 
the only possibilities for n < 7 and the resulting polynomial equations are 

{l + 2w + 4w^ -Aw^){l + 2w-4w'^ -4w^) = 0, for x^^\ and 
(l-w2)(l- W + 16ii;^) = 0, for x'-^l (D.6) 

For i > 2, (i — (a or Cb and (pi = (pa or (pb as the analogue of Case 4. The 
reduction of the singularity conditions to polynomial form is similar to that described 
in the Case 4 analysis and will not be described further. The only possibilities for 
n <7 are 

(1 - 2u; - W + 12w^ + 16w^){l -w + 9w'^ - 2Aw^ + 16w^) 
{l-w + w'^+ 80w^ + 352u;^ + 512u;^ + 256w^) 

(1 + 2u; - 20«;2 - 68u;^ - 32u;^ - 64w^) = 0, for x^^\ and 

(D.7) 

+ 16w'^)(l - 41^2 + 640w7^ - 4096?i;^) 
(16 - 664m;2 + 11273^;'' - 68290w*^ + 141889m;® - 16896m;^° + 65536^^2) 
(1 - 77u;^ + 1898«;^ - 20282u;<' + 107013u;* - 160553w^° + 198A32w^'^ 
-87776u;^^ - 7680u;^^ + 20224u;i^ + 4096^^0) = 0, for x^^^ • 
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Two of the singularity polynomials for x^^^ , l — w^ and 1 — fz;^ + 16w^, are present 
also in the absence of (G^"-') . All other polynomials in (D.6) and (D.7) are new but 
have not been seen in any ODE analysis. Because the integrand singularity surfaces 
X1X2 = I and n^^i = 1 involve the common factor x-i_X2 it is conceivable that the 
and (j)2 integration contours might always be constrained to lie on the same side of the 
two surfaces and never be pinched in between. This would explain the absence of these 
singularities but to resolve such a complicated topological question is not something 
we have attempted. On the other hand, the absence of these singularities is confirmed 
by the ( (j> interchange symmetry discussed in the Case 3 section and the following 
b) results. 

b) van Hove case: sin((/)i) = sm{(j)2) = 0. 

Only the 0i, (j)2 combinations 0, and tt, tt are allowed since the X1X2 = 1 constraint 
requires Ci = ~C2 and this is only possible with cos(Ci) and cos(C2) either both 
l/{2w) — 1 or l/(2u!) + 1. This is analogous to the (f)i = —<j>2 situation of case a) 
which is that the remaining conditions (D.5) and (B.4) just yield the singularities of 
•^("-2) without (G*^")) . The difference between what we have here and in a) is that 
there are no other possibilities: singularities such as (D.6, D.7) are not generated here 
and, by inference from the C <-> interchange symmetry, are not present in a). 

c) The combination (f)i = —(f>2 and C,i = —(,2 with sin((/)i) ^ 0, sin(Ci) ^ 0. 

The (1 = —^2 condition is required by X1X2 = 1. The remaining </>! = —02 by 
sin((^i) • sin(02) — sin{(f)i) ■ sin(C2) = in (D.4). The conditions (D.5) and (B.4) give 
the singularities of x^"^^^ without (G^"^) as in b). 

d) The case sin((j)„) = 0. 

The conditions (D.5) require, for all « > 2, sin((/)i) = and thus ~ or it. 
The corresponding cos{Q) are cos(^*^^^) = 1/(2^) ± 1. The condition X1X2 = 1 
requires Ci = — C2 or cos(Ci) = cos(C2) and, then, from the definition (B.7) that 
cos((/)i) = cos((/)2), or = ±(^2- Note that = (^2 is allowed because (D.4) is no 
longer a constraint on or ^2 once sin(0„) = 0. In the case that (/>i = —02 the phase 
constraint (/>i = mod 27r is satisfied if there are an even number of (pi = w, i > 2, 
terms. In the 0i = 02 case we can set 0i = 4>2 — and thus accommodate an odd 
number of 0i = n, i > 2, terms as well. The final result is that when sin(0„) = 0, 
the van Hove singularities are the Case 3 singularities (B.IO) for x'"~^^ but with odd 
k allowed and when supplemented by Case 5 are: 

Tk{l/2w + l) = T„_2„-fc (l/2w; - 1) , (D.8) 
< k < n — 2m, m = 1, • • • , [ft/2]. 

The situation here has parallels to that in a) and b). First, the odd k polynomials 
in (D.8) are not seen in any ODE analysis and the absence of the odd k singularities 
(D.8) might well have the same explanation as that suggested for the absence of (D.6) 
in a). Second, recall the remarks just preceding (D.4) that the mixed van Hove/pinch 
situation in which sin((^i) ^ 0, sin((^2) 7^ and sin(^j) = 0, i > 2, is not possible 
and hence yields no singularities. But with the interchange ^ ^ 0, wc get exactly the 
situation described by (D.8) for k odd and thus the absence of these singularities is 
also confirmed. 

The four cases a) d) exhaust all possibilities for singularities of the toy x^"^ 

integrals where a single (/i2)^ or (^-12)^ replaces the Fermionic (G^"^) factor in the 
integrand. Furthermore, by performing the Landau analysis of the (n—l) dimensional 
integral analogs of (5) rather that the (2n — 2) dimensional integrals as in [14] we gain 
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a powerful symmetry related tool that enables us to compare results obtained in two 
different ways. Without this tool we can only conclude that the Landau singularity 
conditions (D.6-D.8) are necessary but not sufScient. With this tool we can eliminate 
(D.6-D.8) as singularities of the integrals. This proves the toy x^"^ integrals with 
(/i2)^ or {hi2)'^ have the same singularities as those without these integrand factors. 
We have not investigated systematically what happens when more terms from 
')) are included but such calculations will have many similarities to what has 
been done in the case of just (/i2)^ or (/ii2)^- Furthermore, we have not found any 
obvious candidates for new singularities. Thus we conjecture that the singularities 

of the integrals x'"' are exactly those of the integrals without the Fermionic (G'-"'') 
factor and the complete list of these is the list we have given in Appendix C. 

Appendix E. Singularity exponents at Landau singularities 

To determine the singular behaviour of x'-"'' at Landau singularities requires only a 
local analysis and calculations similar to those already done for the Case 2 circle 
singularities [3, 4]. If only leading exponent information is of interest, the calculations 
simplify dramatically and essentially reduce to a power counting argument. We sketch 
in this appendix the calculations for the cases most analogous to Case 2 and derive 
the formula for the exponent 

p= [{n-mf -3\/2, < m < n, (E.l) 

applicable in many situations. Here n is either the x^""* index in the irreducible Cases 
2-4 or the corresponding index of the x subgraph in the reducible Case 5. The m in 
(E.l) is either the or n — A; as defined in (B.IO) and (B.16). These indices are given 
as subscripts in our polynomial labeling scheme (cf. Appendix C) to allow a direct 
calculation of p using (E.l). For example, p = (3^ + 2^ — 3)/2 = 5 for polynomial ('^53,2) 
or any of the reducible polynomials (^-/V/^53^2) with N = 5 + 2i, i > 0. Note that 
Case 2 polynomials (^n) = (^n„^o) have m = so that (E.l) becomes p = (n^ — 3)/2, 
the formula for the circle singularities derived in [3, 4] as exponents satisfying (E.l), 
which we will refer to as "normal" exponents. We also choose to call the principal 
disc physical exponents p = — 1 at the ferromagnetic point and p = at the anti- 
ferromagnetic point "normal". There are other situations, possibly arising from the 
cancellation in certain expressions after analytic continuation onto different Riemann 
sheets. Because of our lack of understanding in most of these cases we will simply call 
the other exponents "anomalous" . The following is then to be understood as a very 
tentative approach to the singularity exponent problem given that we do not have 
simple general criteria necessary for the "normal" situation to occur. We are in this 
respect ultimately guided by the agreement or disagreement with the numerical work 
of Sections 4 and 5. In this regard note that (E.l) represents the leading exponent 
at a singularity and the derivation below does not exclude exponents greater than p. 
On the other hand, if an exponent less than p is observed a different mechanism is 
required. An alternative explanation is that the observed exponent is associated with 
a solution of the ODE that is not the integral x*^"^ . 

The power counting analysis requires that we distinguish between the two cases 
sin(^i) = i-fi = Oor^ 0. We begin with the latter and determine first the singularities 
of x*^"-* without the Fermionic factor. The singularities in this case come from the zeros 
of 1 — n^^i = 1— exp(i^^j). For the behaviour of the integrand in the vicinity of a zero 
we define e as the deviation of a; — 1 from the Landau singularity value and S(j)i as the 
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deviation of the corresponding phase values. A Taylor expansion of ^ Q yields J^d = 
27rm + E Scb^ {OQ/dcP,) - e/2 ^ V sin(C,) - 1/(4^) E '50? (1 - cos cos </.,)/ sin^ Ci + . . . 
and the stationary condition (B.5) that dC,i/d(j)i is the same for all i = 1... n 
implies that the linear term is proportional to J^^'Pi ^^'^ thus vanishes because of 
the phase constraint in (9). We conclude that the leading terms in 1 — J^x^ are of 
the form + ^ ajJ^? with A and constants. The rescaling Scpi = Stpi^/e puts this 
denominator factor in the x^"^ integral into homogeneous form proportional to e. The 
Jacobian of the transformation from n — 1 variables S(j) to 5ijj yields another overall 
factor of e("~i)/2 so that the final scaling of the integral singularity is e("~3)/2^ ^YiQ,t 
is to say, a singularity power p= (n — 3)/2 which is exactly that observed in all ODE 
analyses of these integrals. 

The inclusion of the Fermionic factor leads to extra powers of e that can be 
determined by counting as follows. In the Case 2 situation treated in [3, 4], there are 
n(n — l)/2 factors of hf^, each of which has an 0(1) denominator and a numerator 
proportional to {6(t>i — Sctyj)"^ or e(5Vi — S'tpj)'^and thus 0(e). The total singularity 
power then becomes p = (n — 3)/2 + n{n — l)/2 = (n^ — 3)/2. For Case 4, to which the 
present argument applies, there are k phases and {n — k) phases (f)}, at a Landau 
singularity. This implies that there are now k{n — k) factors of hfj with numerators 
close to sm^{{4>a — 4'b)/'^) and these are 0(1) and not 0(e). Only the k{k — l)/2 and 
{n— k){n — k— l)/2 factors within each (pa and 0f, set respectively are 0(e). This gives 
the singularity power now as p = (n — 3)/2 + k{k — l)/2 + (n — k){n — k—l)/2 which 
is (E.l) with m = k. 

For the reducible Case 5 a new situation arises in the Fermionic factor. Suppose, 
for example, in our singularity list of {(f), () values there are n — k pairs of type {(pb, Cb) 
but that of the remaining k pairs, one is (— 0a,— Ca) and only A; — 1 are (0a, Ca)- 
Then we know there will be — 1 factors of hfj involving a (—0a, —Co), (0a, Ca) 
combination. The numerator in each of these factors will be close to sin^(0a) and thus 
0(1) and not proportional to (J0i — 50j)^ and 0(e). The denominator on the other 
hand is (1 - XiXjf = [1 - exp(i(Ca + 6Ci) + i(-Ca + SQ))]'' « -{SQ + ^0)' = 0(e). 
In summary, there are fc — 1 Fermionic factor terms each 0(l/e). Furthermore, of 
the original k{k — l)/2 factors involving factors of of type 0a only the remaining 
k{k — l)/2 — (fc — 1) = {k — l){k — 2)/2 are 0(e) as before. Factors involving only 
type 0b or type ±0a and 0t, combinations are also unaffected. Counting all powers 
of e listed above shows that the singularity exponent will be (n — 3)/2 — (fe — 1) + 
(fc - l)(fc - 2)/2 + (n - k)in -k-l)/2=[{k- 2f + (n - kf - 3]/2 which is again 
(E.l) but with m and n — m identified with the irreducible subgraph values k — 2 and 
n — fc. If there is more than one sign reversed pair in the reducible Case 5, then the 
above counting argument can be repeated and leads to the conclusion that (E.l) is 
the general result. 

The arguments when sin(Ci) = i • /i = at the Landau singularity are very 
similar to what is described above and show that (E.l) still applies, in particular for 
Case 3. One technical difference and partial result is worth noting. Each fi is of the 
form •\/(e + ai5(t>i) and thus the appropriate rescaling to make this singular function 
homogeneous in e is 5(j)i — 5ipi ■ e. The singularity power counting in the absence of the 
Fermionic factor then gives an e"^"'^ from the .Jacobian of the transformation, an e^"/^ 
from the n denominator factors fi and another e^^/^ from the denominator 1 — J]^ .x^ 
which reduces, in leading order, to fi- The product of these three e factors is 
£("-3)/2 g^jj^ ^jj^g have as the singularity exponent in the absence of the Fermionic 
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factor the value p= (n — 3)/2 exactly as in the sin(^i) = i • /j ^ situation. 

One set of comparisons of (E.l) with observed exponent values can be found 
in Table 4. While there are many cases of agreement, it is also clear that there are 
"anomalous" exponents that we cannot account for. These include the —2, —7/4, —3/2 
and —5/4 at the ferromagnetic point, and more relevant for the present discussion, the 
powers 2, 5/2 and for the polynomials (^5/^3) and (^5/'^3i,2)- The understanding of 
these cases is left as a challenge for the future. 
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